Le chapitre précédent a présenté les méthodes non paramétriques, comme les k plus proches voisins, qui conservent les données d’entraînement et les consultent au moment de la prédiction. Cette approche est intuitive, mais elle a un coût: les données doivent rester en mémoire, et chaque prédiction requiert de parcourir l’ensemble d’entraînement. Ce chapitre développe une approche différente: plutôt que de garder les données, nous cherchons à les résumer dans un ensemble fixe de paramètres. L’apprentissage devient alors un problème d’optimisation.
Apprentissage supervisé¶
Une ingénieure automobile mesure la distance de freinage d’un véhicule à différentes vitesses. Ses données ressemblent à ceci:
| Vitesse (mph) | Distance (ft) |
|---|---|
| 4 | 2 |
| 7 | 4 |
| 12 | 20 |
| 18 | 56 |
| 24 | 93 |
Elle veut prédire la distance de freinage à 30 mph sans faire l’essai. Pour cela, elle cherche une fonction telle que sur ses observations. Si la fonction capture la relation sous-jacente, elle devrait donner une prédiction raisonnable pour des vitesses non mesurées.
Source
import numpy as np
import matplotlib.pyplot as plt
# Configuration pour des figures haute résolution
%config InlineBackend.figure_format = 'retina'
# Données de freinage (Ezekiel, 1930): vitesse (mph) vs distance d'arrêt (ft)
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)
plt.figure(figsize=(6, 4))
plt.scatter(speed, dist, alpha=0.7, label='Observations')
# Fit quadratic
coeffs = np.polyfit(speed, dist, 2)
speed_grid = np.linspace(0, 30, 100)
dist_pred = np.polyval(coeffs, speed_grid)
plt.plot(speed_grid, dist_pred, 'k--', alpha=0.6, label='Fonction ajustée')
# Prediction at 30 mph
pred_30 = np.polyval(coeffs, 30)
plt.scatter([30], [pred_30], marker='x', s=80, color='C1', zorder=5, label=f'Prédiction à 30 mph: {pred_30:.0f} ft')
plt.xlabel('Vitesse (mph)')
plt.ylabel('Distance de freinage (ft)')
plt.legend()
plt.tight_layout()
Ce processus est l’ajustement de courbe (curve fitting). Nous avons des paires (entrée, sortie), nous ajustons une fonction, et nous utilisons cette fonction pour prédire. L’apprentissage supervisé généralise cette idée: les entrées peuvent être des vecteurs de dimension quelconque, les sorties peuvent être continues ou discrètes, et les fonctions candidates peuvent être bien plus complexes qu’un polynôme.
Formellement, nous disposons d’un jeu de données composé de paires, où chaque est une entrée et est la sortie correspondante. L’objectif est de trouver une fonction qui approxime bien la relation entre entrées et sorties, y compris pour des exemples que nous n’avons pas encore observés.
Dans de nombreuses applications, les entrées sont des vecteurs de caractéristiques. Chaque exemple est un vecteur de dimension , où chaque composante représente une mesure ou un attribut. Pour prédire le prix d’une maison, les entrées pourraient être la superficie, le nombre de chambres et l’âge du bâtiment. Pour filtrer les pourriels, les entrées pourraient être des fréquences de mots. Pour diagnostiquer une maladie, les entrées pourraient être des résultats d’analyses sanguines.
Lorsque la sortie est une valeur continue, nous parlons de régression: pour une sortie scalaire, ou pour une sortie vectorielle. La distance de freinage, le prix d’une maison, la concentration d’un médicament dans le sang sont des exemples de régression.
Lorsque la sortie est une catégorie parmi un ensemble fini, nous parlons de classification. Pour la classification binaire, . Pour la classification multiclasse avec catégories, . Déterminer si un courriel est un pourriel, diagnostiquer une maladie, ou reconnaître un chiffre manuscrit sont des exemples de classification.
Mesurer l’erreur¶
Pour choisir entre deux fonctions candidates, nous avons besoin d’un critère qui quantifie la qualité des prédictions. Une fonction de perte mesure l’écart entre une prédiction et la vraie valeur . Une perte de zéro indique une prédiction parfaite; plus la perte est grande, plus l’erreur est importante.
Pour la régression, nous utilisons généralement la perte quadratique:
Cette perte pénalise les grandes erreurs de manière quadratique. Une erreur de 2 coûte quatre fois plus qu’une erreur de 1.
Reprenons les données de freinage. Supposons que notre fonction prédise 50 ft pour une vitesse où la vraie distance est 56 ft. La perte quadratique est . Si elle prédit 70 ft, la perte est . La perte quadratique pénalise sévèrement les grandes erreurs.
Source
import numpy as np
import matplotlib.pyplot as plt
# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)
coeffs = np.polyfit(speed, dist, 2)
predictions = np.polyval(coeffs, speed)
residuals = dist - predictions
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
# Left: predictions vs observations
ax = axes[0]
ax.scatter(speed, dist, alpha=0.7, label='Observations')
speed_grid = np.linspace(4, 25, 100)
ax.plot(speed_grid, np.polyval(coeffs, speed_grid), 'k--', alpha=0.6, label='Prédictions')
for i in range(len(speed)):
ax.plot([speed[i], speed[i]], [dist[i], predictions[i]], 'C1-', alpha=0.5)
ax.set_xlabel('Vitesse (mph)')
ax.set_ylabel('Distance (ft)')
ax.legend()
ax.set_title('Résidus: écarts entre observations et prédictions')
# Right: histogram of squared residuals
ax = axes[1]
ax.hist(residuals**2, bins=15, edgecolor='black', alpha=0.7)
ax.set_xlabel(r'Perte quadratique $(y - \hat{y})^2$')
ax.set_ylabel('Fréquence')
ax.set_title(f'EQM = {np.mean(residuals**2):.1f}')
plt.tight_layout()
Pour la classification, un choix naturel est la perte 0-1:
Cette perte compte simplement les erreurs: elle vaut 1 pour une mauvaise prédiction, 0 sinon.
Le choix de la fonction de perte dépend du problème. En diagnostic médical, manquer une maladie grave (faux négatif) peut avoir des conséquences bien plus importantes que de prescrire un test supplémentaire à un patient sain (faux positif). Une perte asymétrique refléterait cette différence. En régression, si les grandes erreurs sont particulièrement problématiques, la perte quadratique est appropriée; si nous voulons être robustes aux valeurs aberrantes, la perte absolue est préférable.
Le risque¶
La perte évalue une seule prédiction. Pour évaluer un modèle dans son ensemble, nous voulons mesurer sa performance moyenne sur toutes les données possibles, pas seulement sur les exemples que nous avons observés.
Pourquoi des variables aléatoires?¶
Une question naturelle se pose: si nous ajustons une fonction déterministe à des données, pourquoi avons-nous besoin de variables aléatoires et d’espérances? La fonction obtenue n’est-elle pas simplement une courbe fixe?
La réponse tient en un mot: généralisation. Nous ne nous intéressons pas vraiment à la performance sur les données d’entraînement car ces points sont déjà connus. Ce qui compte, c’est la performance sur des données futures que nous n’avons pas encore observées.
Considérons l’exemple de la distance de freinage. Les 50 mesures dans notre tableau sont un échantillon de toutes les mesures possibles. Si nous retournions sur le terrain et mesurions à nouveau, nous obtiendrions des valeurs légèrement différentes. En effet, le même véhicule à 20 mph ne s’arrête pas exactement à la même distance à chaque essai. Il y a de la variabilité intrinsèque: état de la route, température des freins, réflexes du conducteur.
Cette variabilité est capturée par une distribution de probabilité . Nos 50 points sont des tirages de cette distribution. La question fondamentale devient alors:
Notre modèle prédira-t-il bien sur de nouveaux tirages de cette même distribution?
La fonction elle-même est déterministe une fois entraînée. Mais son évaluation, savoir si elle prédit bien ou mal, dépend de quelles données futures elle rencontrera. Et ces données futures sont incertaines: elles seront tirées de , mais nous ne savons pas lesquelles.
Le risque formalise cette idée: c’est la perte moyenne que subira notre modèle lorsqu’il sera confronté à des données tirées de . C’est une mesure de performance prospective, tournée vers le futur.
Définition formelle¶
Le risque d’une fonction est l’espérance de la perte sur la distribution des données:
Décomposons cette formule étape par étape:
: L’espérance mathématique signifie “moyenne sur tous les exemples possibles”. La notation indique que nous tirons les paires selon la distribution de la nature.
: Pour chaque exemple aléatoire , nous calculons la perte entre la vraie valeur et la prédiction du modèle.
L’intégrale : Cette intégrale calcule une moyenne pondérée. Pour chaque paire possible , nous multiplions la perte par la probabilité que cette paire apparaisse dans la nature, puis nous sommons (intégrons) sur toutes les paires possibles.
Exemple concret¶
Considérons un problème de classification binaire en 2D. Supposons que et . Pour calculer le risque, nous devrions:
Diviser l’espace en une grille fine (par exemple, points)
Pour chaque point de la grille, considérer les deux valeurs possibles de (0 et 1)
Pour chaque combinaison , calculer:
La probabilité que cette combinaison apparaisse
La perte si notre modèle prédit
Faire la somme pondérée:
En pratique, pour un espace continu, cette somme devient une intégrale sur un domaine continu, ce qui est encore plus complexe à calculer.
Visualisons ceci concrètement. La figure suivante montre un problème de classification binaire où chaque classe suit une distribution gaussienne en 2D. Les contours représentent la densité pour chaque classe. La ligne pointillée est la frontière de décision d’un classificateur linéaire. Les régions ombrées indiquent où le classificateur fait des erreurs: la région rouge correspond aux points de classe 0 classés comme classe 1, et la région bleue correspond aux points de classe 1 classés comme classe 0.
Source
import numpy as np
import matplotlib.pyplot as plt
# Paramètres du mélange gaussien (classification binaire 2D)
mu0 = np.array([0.0, 0.0])
mu1 = np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])
def gaussian_pdf(x, mu, cov):
"""PDF d'une gaussienne multivariée."""
d = len(mu)
diff = x - mu
cov_inv = np.linalg.inv(cov)
mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))
# Create grid for visualization
x_range = np.linspace(-3, 5, 200)
y_range = np.linspace(-3, 4, 200)
X_grid, Y_grid = np.meshgrid(x_range, y_range)
pos = np.dstack([X_grid, Y_grid])
# Compute class-conditional densities
p_x_given_0 = gaussian_pdf(pos, mu0, cov)
p_x_given_1 = gaussian_pdf(pos, mu1, cov)
# Joint densities (with equal priors)
prior = 0.5
p_x_y0 = p_x_given_0 * (1 - prior)
p_x_y1 = p_x_given_1 * prior
# Linear decision boundary (Bayes optimal for equal covariances)
# w^T x + b = 0 where w = Sigma^{-1}(mu1 - mu0)
cov_inv = np.linalg.inv(cov)
w = cov_inv @ (mu1 - mu0)
b = -0.5 * (mu1 @ cov_inv @ mu1 - mu0 @ cov_inv @ mu0)
# Decision boundary: w[0]*x + w[1]*y + b = 0 => y = -(w[0]*x + b)/w[1]
x_boundary = np.linspace(-3, 5, 100)
y_boundary = -(w[0] * x_boundary + b) / w[1]
# Classifier prediction: classify as 1 if w^T x + b > 0
predictions = (w[0] * X_grid + w[1] * Y_grid + b) > 0
# Misclassification regions
# Class 0 misclassified as 1: true class is 0, but prediction is 1
misclass_0 = predictions # region where we predict 1
# Class 1 misclassified as 0: true class is 1, but prediction is 0
misclass_1 = ~predictions # region where we predict 0
fig, ax = plt.subplots(figsize=(8, 6))
# Plot class-conditional densities as contours
levels = [0.01, 0.05, 0.1, 0.15]
ax.contour(X_grid, Y_grid, p_x_given_0, levels=levels, colors='C0', alpha=0.7)
ax.contour(X_grid, Y_grid, p_x_given_1, levels=levels, colors='C1', alpha=0.7)
# Shade misclassification regions weighted by probability
# Red: class 0 points incorrectly classified as 1
error_region_0 = np.where(misclass_0, p_x_y0, 0)
# Blue: class 1 points incorrectly classified as 0
error_region_1 = np.where(misclass_1, p_x_y1, 0)
ax.contourf(X_grid, Y_grid, error_region_0, levels=[0.001, 0.01, 0.05, 0.1],
colors=['#ff000010', '#ff000030', '#ff000050'], extend='max')
ax.contourf(X_grid, Y_grid, error_region_1, levels=[0.001, 0.01, 0.05, 0.1],
colors=['#0000ff10', '#0000ff30', '#0000ff50'], extend='max')
# Decision boundary
ax.plot(x_boundary, y_boundary, 'k--', linewidth=2, label='Frontière de décision')
# Class centers
ax.scatter(*mu0, s=100, c='C0', marker='x', linewidths=3, zorder=5, label='Centre classe 0')
ax.scatter(*mu1, s=100, c='C1', marker='x', linewidths=3, zorder=5, label='Centre classe 1')
ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.legend(loc='upper left')
ax.set_aspect('equal')
plt.tight_layout()
Le risque est l’intégrale de la perte sur tout l’espace, pondérée par . Les régions ombrées contribuent au risque: chaque point dans ces régions est mal classé, et sa contribution dépend de la densité de probabilité à cet endroit. Les régions denses proches de la frontière contribuent le plus au risque.
Pourquoi le risque est important¶
Le risque mesure ce que nous obtiendrons en moyenne si nous appliquons à de nouvelles données tirées de la même distribution. Un modèle avec un faible risque fait de bonnes prédictions en général, pas seulement sur les exemples d’entraînement. C’est exactement ce que nous voulons optimiser: un modèle qui performe bien sur des données jamais vues, pas seulement sur celles qu’il a déjà observées.
Cette quantité est ce que nous voulons minimiser. Le problème fondamental est que nous ne connaissons pas la distribution de la nature. Nous n’y avons accès qu’indirectement, via un échantillon fini .
Le risque empirique¶
Puisque le risque est inaccessible, nous l’approximons par une moyenne sur les données disponibles. Le risque empirique est:
Cette quantité est calculable: c’est la moyenne des pertes sur l’échantillon d’entraînement. Pour la perte 0-1, le risque empirique est le taux d’erreur sur les données d’entraînement. Pour la perte quadratique, c’est l’erreur quadratique moyenne.
Pourquoi le risque est-il inaccessible?¶
La nécessité d’utiliser le risque empirique découle de deux obstacles fondamentaux, l’un conceptuel et l’autre computationnel.
Obstacle 1: La distribution est inconnue¶
La nature possède une distribution qui génère les données, mais nous ne la connaissons pas. Nous n’observons qu’un échantillon fini tiré de cette distribution.
L’ensemble est une variable aléatoire: si nous répétions l’expérience de collecte de données, nous obtiendrions un échantillon différent. Cette perspective, adoptée notamment dans Murphy (2022), rappelle que nos conclusions dépendent de l’échantillon particulier que nous avons observé. C’est comme si nous regardions quelques gouttes d’eau d’un océan: nous pouvons analyser ces gouttes, mais un autre prélèvement donnerait des gouttes différentes.
Même si nous tentions d’estimer à partir des données (par exemple, via des techniques d’estimation de densité comme les mélanges de gaussiennes ou les estimateurs à noyau), nous n’obtiendrions qu’une approximation de la vraie distribution. Cette approximation serait elle-même imparfaite et dépendrait de nos hypothèses sur la forme de la distribution.
La figure suivante illustre ce problème. À gauche, la vraie distribution que la nature utilise pour générer les données (que nous ne connaissons pas). À droite, un échantillon de points tirés de cette distribution (ce que nous observons).
Source
import numpy as np
import matplotlib.pyplot as plt
# Paramètres du mélange gaussien
mu0, mu1 = np.array([0.0, 0.0]), np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])
def gaussian_pdf(x, mu, cov):
d = len(mu)
diff = x - mu
cov_inv = np.linalg.inv(cov)
mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))
# Générer un échantillon
rng = np.random.default_rng(42)
n = 50
X0 = rng.multivariate_normal(mu0, cov, n // 2)
X1 = rng.multivariate_normal(mu1, cov, n // 2)
X = np.vstack([X0, X1])
y = np.concatenate([np.zeros(n // 2), np.ones(n // 2)])
# Grille pour visualisation
x_range = np.linspace(-3, 5, 150)
y_range = np.linspace(-3, 4, 150)
X_grid, Y_grid = np.meshgrid(x_range, y_range)
pos = np.dstack([X_grid, Y_grid])
p_x_given_0 = gaussian_pdf(pos, mu0, cov)
p_x_given_1 = gaussian_pdf(pos, mu1, cov)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Left: True distribution (what nature knows)
ax = axes[0]
ax.contourf(X_grid, Y_grid, p_x_given_0, levels=15, cmap='Blues', alpha=0.6)
ax.contourf(X_grid, Y_grid, p_x_given_1, levels=15, cmap='Oranges', alpha=0.6)
ax.contour(X_grid, Y_grid, p_x_given_0, levels=5, colors='C0', alpha=0.8)
ax.contour(X_grid, Y_grid, p_x_given_1, levels=5, colors='C1', alpha=0.8)
ax.scatter(*mu0, s=100, c='C0', marker='x', linewidths=3, zorder=5)
ax.scatter(*mu1, s=100, c='C1', marker='x', linewidths=3, zorder=5)
ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('Distribution vraie $p(x, y)$\n(inconnue)')
ax.set_aspect('equal')
# Right: Finite sample (what we observe)
ax = axes[1]
ax.scatter(X[y == 0, 0], X[y == 0, 1], c='C0', alpha=0.7, s=50, label='Classe 0')
ax.scatter(X[y == 1, 0], X[y == 1, 1], c='C1', alpha=0.7, s=50, label='Classe 1')
ax.set_xlim(-3, 5)
ax.set_ylim(-3, 4)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title(f'Échantillon observé $\\mathcal{{D}}$\n($N = {len(X)}$ points)')
ax.legend()
ax.set_aspect('equal')
plt.tight_layout()
Nous ne voyons que les points à droite. La structure continue à gauche, incluant les contours, les densités, ainsi que les régions de haute et basse probabilité, nous est cachée. C’est à partir de ces quelques points que nous devons estimer la performance de notre modèle.
Obstacle 2: L’intégration est computationnellement intractable¶
Supposons, par un miracle, que nous connaissions exactement . Pourrions-nous alors calculer le risque ?
La réponse est généralement non, pour plusieurs raisons:
Pour les espaces continus: L’intégrale est une intégrale de grande dimension. Si avec grand (par exemple, pour des images ou pour des données textuelles), nous devons intégrer sur un espace de dimension .
Pour vous rappeler l’idée, en calcul on approche une intégrale en 1D par une somme: on découpe l’intervalle en petites tranches et on additionne des aires de rectangles ou de trapèzes. Par exemple, sur :
Cette idée générale, qui consiste à remplacer une intégrale par une somme pondérée de valeurs de évaluées à des points , s’appelle l’intégration numérique (ou quadrature).
Le problème en apprentissage est que notre intégrale n’est pas en 1D. Si on applique le même raisonnement en dimension en mettant, disons, points par dimension, on obtient une grille de taille (et ici peut être très grand). Le nombre de points à évaluer explose donc exponentiellement avec . C’est exactement la malédiction de la dimensionnalité.
Pour les espaces discrets: Si et sont discrets mais prennent de nombreuses valeurs, la somme peut avoir un nombre exponentiel de termes. Par exemple, si est un vecteur binaire de dimension , il y a valeurs possibles pour . Pour , cela fait déjà termes à sommer, ce qui est computationnellement impossible.
Intégration de Monte Carlo: On pourrait penser utiliser l’intégration de Monte Carlo: tirer des échantillons selon et estimer l’intégrale par la moyenne empirique. Mais pour obtenir une estimation précise du risque, nous aurions besoin d’un très grand nombre d’échantillons (potentiellement infini pour une précision parfaite). De plus, cela nécessiterait de pouvoir échantillonner efficacement depuis , ce qui est lui-même un problème difficile si la distribution est complexe.
La figure suivante illustre la malédiction de la dimensionnalité. Avec seulement 10 points par dimension pour une quadrature numérique, le nombre total de points d’évaluation explose rapidement.
Source
import numpy as np
import matplotlib.pyplot as plt
# Number of grid points per dimension
points_per_dim = 10
# Dimensions to consider
dimensions = np.array([1, 2, 3, 5, 10, 20, 50, 100])
# Total grid points = points_per_dim^d
total_points = points_per_dim ** dimensions.astype(float)
fig, ax = plt.subplots(figsize=(8, 5))
bars = ax.bar(range(len(dimensions)), total_points, color='steelblue', edgecolor='black')
# Add reference lines
ax.axhline(y=1e9, color='C1', linestyle='--', alpha=0.7, label='1 milliard (limite pratique)')
ax.axhline(y=1e80, color='C3', linestyle=':', alpha=0.7, label='$10^{80}$ (atomes dans l\'univers)')
ax.set_yscale('log')
ax.set_xticks(range(len(dimensions)))
ax.set_xticklabels([f'd={d}' for d in dimensions])
ax.set_xlabel('Dimension de l\'espace des entrées')
ax.set_ylabel('Nombre de points de grille')
ax.set_title(f'Points nécessaires pour l\'intégration numérique\n({points_per_dim} points par dimension)')
ax.legend(loc='upper left')
# Annotate a few bars
for i, (d, n) in enumerate(zip(dimensions, total_points)):
if d <= 5:
ax.annotate(f'$10^{{{d}}}$', (i, n), ha='center', va='bottom', fontsize=9)
elif d == 10:
ax.annotate(f'$10^{{{10}}}$', (i, n), ha='center', va='bottom', fontsize=9)
elif d == 100:
ax.annotate(f'$10^{{{100}}}$', (i, n), ha='center', va='bottom', fontsize=9)
ax.set_ylim(1, 1e105)
plt.tight_layout()
En dimension 10, il faut déjà 1010 points, soit dix milliards. En dimension 100, il en faut 10100, un nombre qui dépasse le nombre d’atomes dans l’univers observable. L’intégration numérique directe est donc impossible pour les problèmes de haute dimension, même si nous connaissions exactement.
Le risque empirique comme seule option pratique¶
Face à ces obstacles, le risque empirique est notre seule option calculable. Mais il y a une bonne nouvelle: le risque empirique est une forme d’intégration de Monte Carlo, et Monte Carlo a une propriété remarquable.
| Méthode | Complexité | Exigence |
|---|---|---|
| Quadrature (règles trapézoïdales, etc.) | Connaître exactement | |
| Monte Carlo | Avoir des échantillons de |
La complexité de Monte Carlo est indépendante de la dimension . Elle ne dépend que du nombre d’échantillons . C’est cette propriété qui rend l’apprentissage possible en haute dimension. De plus, nous n’avons pas besoin de connaître la valeur numérique de : nous avons seulement besoin de pouvoir tirer des échantillons de cette distribution. C’est exactement ce que nos données d’entraînement nous fournissent.
Le risque empirique remplace l’intégrale sur la distribution inconnue par une moyenne sur l’échantillon fini que nous possédons:
Cette formule est directe à évaluer: nous parcourons nos exemples d’entraînement, calculons la perte pour chacun, et faisons la moyenne.
Reprenons les données de freinage. Divisons-les en deux parties: les mesures à vitesses faibles (4-19 mph) pour l’entraînement, et les mesures à vitesses élevées (20-25 mph) pour le test. Le risque empirique sur l’ensemble d’entraînement mesure la qualité de l’ajustement. Le risque empirique sur l’ensemble de test estime la performance sur des vitesses non vues pendant l’entraînement.
Source
import numpy as np
import matplotlib.pyplot as plt
# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)
# Split: train on low speeds, test on high speeds
train_mask = speed < 20
test_mask = speed >= 20
speed_train, dist_train = speed[train_mask], dist[train_mask]
speed_test, dist_test = speed[test_mask], dist[test_mask]
# Fit on training data
coeffs = np.polyfit(speed_train, dist_train, 2)
# Compute MSE on train and test
pred_train = np.polyval(coeffs, speed_train)
pred_test = np.polyval(coeffs, speed_test)
mse_train = np.mean((dist_train - pred_train)**2)
mse_test = np.mean((dist_test - pred_test)**2)
plt.figure(figsize=(7, 4))
plt.scatter(speed_train, dist_train, alpha=0.7, label=f'Entraînement (EQM={mse_train:.1f})')
plt.scatter(speed_test, dist_test, alpha=0.7, marker='s', label=f'Test (EQM={mse_test:.1f})')
speed_grid = np.linspace(4, 28, 100)
plt.plot(speed_grid, np.polyval(coeffs, speed_grid), 'k--', alpha=0.6, label='Fonction ajustée')
plt.axvline(x=20, color='gray', linestyle=':', alpha=0.5)
plt.xlabel('Vitesse (mph)')
plt.ylabel('Distance (ft)')
plt.legend()
plt.tight_layout()
Dans cet exemple, l’EQM (MSE, erreur quadratique moyenne) sur l’ensemble de test est plus élevé que sur l’ensemble d’entraînement. Cet écart est typique: la fonction a été optimisée pour les données d’entraînement, pas pour les données de test.
Sous l’hypothèse que les exemples sont tirés indépendamment et identiquement distribués (i.i.d.) selon , le risque empirique est un estimateur non biaisé du vrai risque: . Cela signifie qu’en moyenne, sur tous les échantillons possibles, le risque empirique est égal au vrai risque.
Par la loi des grands nombres, lorsque , le risque empirique converge vers le vrai risque (presque sûrement). Avec suffisamment de données, si l’échantillon est représentatif de la distribution, le risque empirique devrait être proche du risque.
La figure suivante illustre cette convergence. Nous utilisons le problème de classification gaussienne pour lequel nous pouvons calculer le vrai risque analytiquement. Chaque courbe montre l’évolution du risque empirique pour un échantillon de taille croissante. Toutes les courbes convergent vers le vrai risque (ligne pointillée), mais avec des fluctuations qui diminuent à mesure que augmente.
Source
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# Paramètres du mélange gaussien
mu0, mu1 = np.array([0.0, 0.0]), np.array([2.0, 1.0])
cov = np.array([[1.0, 0.3], [0.3, 1.0]])
def gaussian_pdf(x, mu, cov):
d = len(mu)
diff = x - mu
cov_inv = np.linalg.inv(cov)
mahal = np.einsum('...i,ij,...j->...', diff, cov_inv, diff)
return np.exp(-0.5 * mahal) / np.sqrt((2 * np.pi) ** d * np.linalg.det(cov))
# Compute Bayes-optimal classifier error rate (true risk)
# For Gaussian classes with equal covariance, the Bayes error is:
# P(error) = Phi(-d/2) where d is the Mahalanobis distance between means
cov_inv = np.linalg.inv(cov)
d_squared = (mu1 - mu0) @ cov_inv @ (mu1 - mu0)
d = np.sqrt(d_squared)
true_risk = norm.cdf(-d / 2)
# Simulate empirical risk for different sample sizes
sample_sizes = np.arange(10, 1001, 10)
n_runs = 20
fig, ax = plt.subplots(figsize=(9, 5))
# Store all runs for confidence band
all_risks = np.zeros((n_runs, len(sample_sizes)))
for run in range(n_runs):
empirical_risks = []
# Generate a large dataset and compute cumulative empirical risk
rng = np.random.default_rng(run)
n_total = 1000
X0 = rng.multivariate_normal(mu0, cov, n_total // 2)
X1 = rng.multivariate_normal(mu1, cov, n_total // 2)
X = np.vstack([X0, X1])
y = np.concatenate([np.zeros(n_total // 2), np.ones(n_total // 2)])
perm = rng.permutation(n_total)
X, y = X[perm], y[perm]
# Bayes-optimal classifier: predict 1 if w^T x + b > 0
w = cov_inv @ (mu1 - mu0)
b = -0.5 * (mu1 @ cov_inv @ mu1 - mu0 @ cov_inv @ mu0)
for n in sample_sizes:
X_n, y_n = X[:n], y[:n]
predictions = (X_n @ w + b > 0).astype(float)
emp_risk = np.mean(predictions != y_n)
empirical_risks.append(emp_risk)
all_risks[run] = empirical_risks
ax.plot(sample_sizes, empirical_risks, 'C0-', alpha=0.15, linewidth=0.8)
# Mean and confidence bands
mean_risk = np.mean(all_risks, axis=0)
std_risk = np.std(all_risks, axis=0)
ax.fill_between(sample_sizes, mean_risk - 2*std_risk, mean_risk + 2*std_risk,
alpha=0.3, color='C0', label='Intervalle ± 2 écarts-types')
ax.plot(sample_sizes, mean_risk, 'C0-', linewidth=2, label='Moyenne empirique')
# True risk
ax.axhline(y=true_risk, color='C3', linestyle='--', linewidth=2,
label=f'Vrai risque = {true_risk:.3f}')
ax.set_xlabel('Taille de l\'échantillon $N$')
ax.set_ylabel('Risque empirique (taux d\'erreur)')
ax.set_title('Convergence du risque empirique vers le vrai risque')
ax.legend(loc='upper right')
ax.set_xlim(0, 1000)
ax.set_ylim(0, 0.35)
plt.tight_layout()
Avec , le risque empirique peut facilement varier de 0.10 à 0.25 selon l’échantillon. Avec , la variabilité est beaucoup plus faible. C’est la loi des grands nombres en action: plus l’échantillon est grand, plus l’estimation est précise.
Le compromis fondamental¶
Cette situation crée un compromis fondamental en apprentissage automatique:
Ce que nous voulons minimiser: Le risque , qui mesure la performance sur toutes les données possibles
Ce que nous pouvons minimiser: Le risque empirique , qui mesure la performance sur nos données d’entraînement
L’écart entre ces deux quantités est au cœur de l’apprentissage automatique. Un modèle peut avoir un risque empirique très faible (il performe bien sur les données d’entraînement) tout en ayant un risque élevé (il performe mal sur de nouvelles données). C’est le problème du surapprentissage, que nous explorerons plus en détail dans le chapitre sur la généralisation.
La question de savoir quand et à quelle vitesse l’approximation du risque par le risque empirique est fiable relève de la théorie de la généralisation, que nous aborderons au chapitre suivant.
Minimisation du risque empirique¶
Nous avons maintenant les éléments pour formuler l’apprentissage comme un problème d’optimisation. Nous cherchons la fonction dans une classe qui minimise le risque:
Puisque le risque est inaccessible, nous le remplaçons par le risque empirique:
Ce principe est la minimisation du risque empirique (MRE): choisir la fonction qui fait le moins d’erreurs sur les données d’entraînement, en espérant que cette performance se transfère aux nouvelles données.
La classe est notre classe d’hypothèses. Elle représente l’ensemble des fonctions que nous sommes prêts à considérer. Le choix de encode nos hypothèses sur la forme de la relation entre entrées et sorties.
Un premier exemple: les modèles linéaires¶
Pour rendre ces concepts concrets, considérons la classe la plus simple: les modèles linéaires. Un modèle linéaire suppose que la sortie est une combinaison linéaire des entrées:
où est le vecteur d’entrée augmenté d’un 1 pour le biais (), et est le vecteur de paramètres contenant le biais et les poids .
Cette forme est restrictive: elle suppose que la relation entre entrées et sorties est linéaire. Pour les données de freinage, cela signifierait que la distance est proportionnelle à la vitesse, ce qui n’est pas le cas (la relation est plutôt quadratique). Néanmoins, les modèles linéaires sont utiles comme point de départ, et nous verrons comment les étendre pour capturer des relations non linéaires.
Avec cette classe fixée, l’apprentissage consiste à trouver les paramètres qui minimisent le risque empirique. Pour la perte quadratique, cela revient à minimiser la somme des carrés des résidus.
Mais quand le minimiseur du risque empirique a-t-il un faible risque? Si minimise et minimise , nous voulons que soit proche de . La réponse dépend de la taille de l’échantillon , de la complexité de la classe , et de propriétés de la distribution .
Résoudre le problème d’optimisation¶
Quand nous écrivons , nous cherchons les paramètres qui rendent la fonction objectif aussi petite que possible. Mais comment trouver ces paramètres en pratique?
La réponse dépend de la forme du problème:
Pour certains problèmes, comme la régression linéaire avec perte quadratique, nous pouvons dériver une solution analytique en posant le gradient égal à zéro et en résolvant le système d’équations résultant.
Pour d’autres, nous devons recourir à des algorithmes itératifs (comme la descente de gradient) ou à des solveurs spécialisés (comme la programmation quadratique pour les SVM).
Exemple: solution analytique pour la régression linéaire (MCO)¶
Pour illustrer les solutions analytiques, considérons la régression linéaire avec perte quadratique. L’objectif est de minimiser la somme des carrés des résidus:
où est la matrice des entrées (avec une colonne de 1 pour le biais) et est le vecteur des sorties.
En développant et en calculant le gradient:
En posant le gradient égal à zéro, nous obtenons les équations normales:
Si la matrice est inversible, la solution unique est:
Cette solution porte le nom d’estimateur des moindres carrés ordinaires (MCO, ou ordinary least squares, OLS). Elle peut être calculée directement sans itération, ce qui en fait un exemple classique de solution analytique.
Le chapitre sur l’optimisation présentera ces méthodes en détail. Pour l’instant, nous nous concentrons sur la formulation des problèmes d’apprentissage comme problèmes d’optimisation, en gardant à l’esprit que des outils existent pour les résoudre.
Généralisation¶
La différence entre le risque et le risque empirique est l’écart de généralisation:
Un modèle qui minimise le risque empirique peut avoir un risque élevé si cet écart est grand. Ce phénomène est le surapprentissage: le modèle s’ajuste aux particularités de l’échantillon d’entraînement, y compris le bruit, plutôt qu’aux régularités sous-jacentes. L’erreur d’entraînement est faible, mais l’erreur sur de nouvelles données est élevée.
À l’inverse, un modèle trop simple peut avoir un risque empirique et un risque tous deux élevés. C’est le sous-apprentissage: le modèle n’a pas la capacité de capturer la structure des données.
Extrapolation¶
Un cas particulier de mauvaise généralisation est l’extrapolation: prédire pour des entrées en dehors de la plage des données d’entraînement. Même un modèle bien ajusté peut échouer spectaculairement lorsqu’on lui demande de prédire au-delà de ce qu’il a vu.
Considérons des essais en soufflerie pour mesurer la portance d’une aile à différentes vitesses. Les tests sont effectués entre 20 et 60 m/s. L’ingénieur veut prédire la portance à 100 m/s.
Source
import numpy as np
import matplotlib.pyplot as plt
# Données de portance aérodynamique (simulées)
np.random.seed(42)
rho, S, C_L = 1.225, 20.0, 0.5
v_train = np.linspace(20, 60, 8)
L_true_train = 0.5 * rho * v_train**2 * S * C_L
L_train = L_true_train + np.random.normal(0, 400, len(v_train))
coeffs_2 = np.polyfit(v_train, L_train, 2)
coeffs_5 = np.polyfit(v_train, L_train, 5)
v_extrap = np.linspace(15, 110, 200)
L_true_extrap = 0.5 * rho * v_extrap**2 * S * C_L
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
for ax, coeffs, deg in zip(axes, [coeffs_2, coeffs_5], [2, 5]):
ax.scatter(v_train, L_train, s=50, zorder=5, label='Observations')
ax.plot(v_extrap, L_true_extrap, 'g-', alpha=0.5, label='Vraie relation')
L_pred = np.polyval(coeffs, v_extrap)
ax.plot(v_extrap, L_pred, 'k--', label=f'Polynôme degré {deg}')
ax.axvline(60, color='gray', linestyle=':', alpha=0.5)
ax.axvspan(60, 110, alpha=0.1, color='red')
ax.set_xlabel('Vitesse (m/s)')
ax.set_ylabel('Portance (N)')
ax.set_title(f'Degré {deg}')
ax.legend(loc='upper left')
ax.set_ylim(-5000, 80000)
ax.text(85, 5000, 'Extrapolation', ha='center', fontsize=10, color='red', alpha=0.7)
plt.tight_layout()
Le polynôme de degré 2 (qui correspond au vrai modèle physique ) extrapole correctement. Le polynôme de degré 5, bien qu’il ajuste aussi bien les données d’entraînement, diverge complètement en dehors de la plage observée.
Régularisation¶
Une manière de contrôler le surapprentissage consiste à pénaliser la complexité du modèle directement dans la fonction objectif. Le risque empirique régularisé est:
où mesure la complexité du modèle et contrôle l’intensité de la pénalisation. Un choix courant est la régularisation (ou weight decay):
Cette pénalisation pousse les paramètres vers zéro, ce qui a pour effet de lisser la fonction apprise. En régression linéaire, l’ajout de cette pénalité donne la régression ridge:
Illustrons l’effet de la régularisation sur le même problème de régression polynomiale. Avec un polynôme de degré 15 et différentes valeurs de :
Source
import numpy as np
import matplotlib.pyplot as plt
# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)
# Train/test split
np.random.seed(42)
indices = np.random.permutation(len(speed))
train_idx, test_idx = indices[:35], indices[35:]
speed_train, dist_train = speed[train_idx], dist[train_idx]
speed_test, dist_test = speed[test_idx], dist[test_idx]
# Build polynomial features (degree 15)
degree = 15
def poly_features(x, deg):
return np.vstack([x**i for i in range(deg+1)]).T
X_train = poly_features(speed_train, degree)
X_test = poly_features(speed_test, degree)
# Ridge regression for different lambda values
lambdas = [0, 1e-6, 1e-3, 1]
fig, axes = plt.subplots(2, 2, figsize=(10, 8))
for ax, lam in zip(axes.flat, lambdas):
# Ridge solution: (X^T X + lambda I)^{-1} X^T y
I = np.eye(X_train.shape[1])
I[0, 0] = 0 # Don't regularize bias
w = np.linalg.solve(X_train.T @ X_train + lam * I, X_train.T @ dist_train)
# Predictions
pred_train = X_train @ w
pred_test = X_test @ w
mse_train = np.mean((dist_train - pred_train)**2)
mse_test = np.mean((dist_test - pred_test)**2)
# Plot
ax.scatter(speed_train, dist_train, alpha=0.6, s=30, label='Entraînement')
ax.scatter(speed_test, dist_test, alpha=0.6, s=30, marker='s', label='Test')
speed_grid = np.linspace(3, 26, 200)
X_grid = poly_features(speed_grid, degree)
pred_grid = X_grid @ w
pred_grid = np.clip(pred_grid, -50, 200)
ax.plot(speed_grid, pred_grid, 'k-', alpha=0.7)
ax.set_xlim(3, 26)
ax.set_ylim(-20, 150)
ax.set_xlabel('Vitesse (mph)')
ax.set_ylabel('Distance (ft)')
ax.set_title(f'$\\lambda$ = {lam}: Entr. EQM={mse_train:.1f}, Test EQM={mse_test:.1f}')
if lam == 0:
ax.legend()
plt.tight_layout()
Sans régularisation (), le polynôme de degré 15 oscille fortement. Avec une régularisation modérée (), les oscillations sont atténuées et l’erreur de test diminue. Avec une régularisation trop forte (), le modèle devient trop contraint et sous-apprend.
Solution analytique de la régression ridge¶
Comme pour les moindres carrés ordinaires, la régression ridge admet une solution analytique. L’objectif régularisé est:
En développant et en calculant le gradient:
En posant le gradient égal à zéro, nous obtenons les équations normales régularisées:
La solution est:
Comparons avec la solution MCO: . La seule différence est l’ajout du terme à la matrice .
Solution via décomposition en valeurs singulières (SVD)¶
Cette section présente une autre façon d’exprimer les solutions MCO et Ridge, en utilisant la décomposition en valeurs singulières (SVD). Si vous n’avez jamais rencontré la SVD, ne vous inquiétez pas: nous allons l’introduire progressivement. Cette approche n’est pas strictement nécessaire pour comprendre la régression ridge, mais elle offre une interprétation géométrique très éclairante qui révèle pourquoi la régularisation fonctionne.
Qu’est-ce que la SVD? Si vous avez déjà rencontré la décomposition en valeurs propres, la SVD en est une généralisation. Pour une matrice carrée symétrique , la décomposition en valeurs propres s’écrit , où contient les vecteurs propres et les valeurs propres. La SVD généralise cette idée à n’importe quelle matrice, même rectangulaire.
Pour une matrice de données, la SVD la réécrit comme le produit de trois matrices:
Lien avec la décomposition en valeurs propres: Les colonnes de sont les vecteurs propres de , et les valeurs singulières sont les racines carrées des valeurs propres de . Autrement dit, si est la décomposition en valeurs propres de , alors où sont les valeurs propres. De même, les colonnes de sont les vecteurs propres de .
Cette connexion est utile car apparaît naturellement dans la régression (c’est la matrice que nous inversons pour MCO). Les valeurs singulières nous renseignent donc directement sur le “conditionnement” de cette matrice: si certaines valeurs singulières sont très petites, alors est proche d’être singulière (non inversible).
Interprétation géométrique:
contient les directions principales dans l’espace des caractéristiques (les colonnes sont orthonormales). Ces directions correspondent aux axes le long desquels la matrice transforme les vecteurs de manière la plus efficace.
est une matrice diagonale contenant les valeurs singulières , ordonnées du plus grand au plus petit. Chaque valeur singulière mesure l’amplitude de la transformation le long de la direction . Une grande valeur singulière signifie que la transformation est forte dans cette direction; une petite valeur singulière signifie que la transformation est faible.
contient les directions correspondantes dans l’espace des observations (les colonnes sont orthonormales). Chaque indique comment les observations se projettent sur la direction principale .
Solution MCO via SVD: En utilisant cette décomposition, la solution MCO peut s’écrire:
Cette formule décompose la solution en une somme de contributions le long de chaque direction principale . Le terme mesure combien la sortie s’aligne avec la direction , divisé par l’amplitude de cette direction. Notez que diviser par une petite valeur singulière peut amplifier le bruit, ce qui explique pourquoi MCO peut être instable quand certaines directions ont de petites valeurs singulières.
Solution Ridge via SVD: Pour Ridge, la solution devient:
La différence avec MCO est le facteur de rétrécissement qui multiplie chaque terme. Ce facteur est toujours inférieur à 1, ce qui “rétrécit” chaque composante vers zéro. L’effet clé est que ce rétrécissement est différencié: les directions avec de petites valeurs singulières sont rétrécies plus fortement que celles avec de grandes valeurs singulières.
Interprétation géométrique: Les directions principales définissent les axes d’une ellipse de confiance dans l’espace des coefficients. Pour MCO, les longueurs des demi-axes sont proportionnelles à : plus une direction a une petite valeur singulière , plus la variance d’estimation est grande (nous verrons cette distinction cruciale plus loin). Avec Ridge, ces longueurs deviennent proportionnelles à : l’ellipse se rétrécit globalement, et plus rapidement le long des directions associées aux petites valeurs singulières. C’est exactement ce que nous voulons: réduire la variance d’estimation là où elle est la plus grande, c’est-à-dire dans les directions où les données sont peu dispersées.
Avantages numériques: Au-delà de l’interprétation, la SVD offre aussi des avantages pratiques. Elle est plus stable numériquement que l’inversion directe de , surtout quand cette matrice est mal conditionnée (c’est-à-dire quand certaines valeurs singulières sont très petites). Les algorithmes SVD gèrent mieux ces cas délicats.
Visualisation: ellipse des données et vecteurs singuliers¶
Pour rendre ces concepts concrets, visualisons ce que la SVD capture sur un nuage de données 2D. Générons des points suivant une distribution gaussienne avec une covariance non triviale (les deux variables sont corrélées).
Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
np.random.seed(42)
# Générer des données gaussiennes corrélées
n_points = 200
mean = [0, 0]
cov = [[2.0, 1.2], [1.2, 1.0]] # Covariance non diagonale
X = np.random.multivariate_normal(mean, cov, n_points)
# Centrer les données
X_centered = X - X.mean(axis=0)
# SVD de la matrice de données centrée
U, d, Vt = np.linalg.svd(X_centered, full_matrices=False)
V = Vt.T
# Les valeurs singulières sont liées aux écarts-types: d_j / sqrt(N-1)
# Pour l'ellipse, nous utilisons les écarts-types dans chaque direction
std_1 = d[0] / np.sqrt(n_points - 1)
std_2 = d[1] / np.sqrt(n_points - 1)
# Créer la figure
fig, ax = plt.subplots(figsize=(8, 6))
# Tracer les points
ax.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.5, s=20, c='tab:blue', label='Données')
# Tracer les vecteurs singuliers (directions principales)
origin = [0, 0]
scale = 2 # Facteur d'échelle pour la visualisation
# Premier vecteur singulier (direction de plus grande variance)
ax.annotate('', xy=V[:, 0] * std_1 * scale, xytext=origin,
arrowprops=dict(arrowstyle='->', color='tab:red', lw=2.5))
ax.annotate('', xy=-V[:, 0] * std_1 * scale, xytext=origin,
arrowprops=dict(arrowstyle='->', color='tab:red', lw=2.5))
# Deuxième vecteur singulier (direction de plus petite variance)
ax.annotate('', xy=V[:, 1] * std_2 * scale, xytext=origin,
arrowprops=dict(arrowstyle='->', color='tab:orange', lw=2.5))
ax.annotate('', xy=-V[:, 1] * std_2 * scale, xytext=origin,
arrowprops=dict(arrowstyle='->', color='tab:orange', lw=2.5))
# Ellipse de confiance (2 écarts-types)
angle = np.degrees(np.arctan2(V[1, 0], V[0, 0]))
ellipse = Ellipse(xy=(0, 0), width=4*std_1, height=4*std_2, angle=angle,
fill=False, edgecolor='gray', linestyle='--', linewidth=1.5)
ax.add_patch(ellipse)
# Annotations
ax.text(V[0, 0] * std_1 * scale * 1.15, V[1, 0] * std_1 * scale * 1.15,
f'$\\mathbf{{v}}_1$ ($d_1 = {d[0]:.1f}$)', fontsize=11, color='tab:red')
ax.text(V[0, 1] * std_2 * scale * 1.3, V[1, 1] * std_2 * scale * 1.3,
f'$\\mathbf{{v}}_2$ ($d_2 = {d[1]:.1f}$)', fontsize=11, color='tab:orange')
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('Nuage gaussien et directions principales (SVD)')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.axhline(0, color='gray', linewidth=0.5)
ax.axvline(0, color='gray', linewidth=0.5)
ax.set_xlim(-4, 4)
ax.set_ylim(-3, 3)
plt.tight_layout()
La figure montre un nuage de 200 points tirés d’une gaussienne 2D. Les flèches représentent les vecteurs singuliers et :
Le vecteur (rouge) pointe dans la direction de plus grande variance. La valeur singulière mesure l’amplitude de la dispersion dans cette direction.
Le vecteur (orange) pointe dans la direction de plus petite variance, perpendiculaire à . La valeur singulière est plus petite.
L’ellipse en pointillés représente la région contenant environ 95% des données si elles suivent exactement la distribution gaussienne. Ses axes coïncident avec les vecteurs singuliers, et les longueurs des demi-axes sont proportionnelles aux valeurs singulières.
Cette visualisation illustre pourquoi la SVD est si utile: elle identifie automatiquement les axes naturels des données. Dans le contexte de la régression, si les caractéristiques forment un nuage allongé (valeurs singulières très différentes), alors certaines directions contiennent beaucoup d’information (grandes valeurs singulières) tandis que d’autres en contiennent peu (petites valeurs singulières).
Deux variances: données vs estimation¶
Avant d’aller plus loin, clarifions une source fréquente de confusion. Le mot variance apparaît dans deux contextes très différents lorsqu’on parle de SVD et de régularisation:
Variance des données (dispersion): mesure l’étalement des données le long d’une direction . Elle est proportionnelle à . Une grande valeur singulière signifie que les données sont très dispersées dans cette direction.
Variance d’estimation (incertitude): mesure l’incertitude sur notre estimé du paramètre correspondant à la direction . Elle est proportionnelle à . Une petite valeur singulière signifie que notre estimé est très incertain.
Ces deux variances sont inversement reliées:
| Valeur singulière | Variance des données | Variance d’estimation | Interprétation |
|---|---|---|---|
| Grande | Élevée (données étalées) | Faible (estimé précis) | Beaucoup d’information |
| Petite | Faible (données concentrées) | Élevée (estimé incertain) | Peu d’information |
Intuition: Imaginez estimer une pente à partir de données. Si les points sont très étalés horizontalement (grande variance des données en ), la pente est facile à déterminer avec précision (faible variance d’estimation). Si les points sont tous regroupés (petite variance des données), la pente est très incertaine (grande variance d’estimation).
C’est cette relation inverse qui explique le comportement de Ridge:
Ridge rétrécit les directions où est petit (faible variance des données)
Ce sont précisément les directions où la variance d’estimation est grande
En rétrécissant ces directions, Ridge réduit la variance d’estimation au prix d’un biais
Ainsi, quand nous disons que « Ridge contrôle la variance », nous parlons de la variance d’estimation des paramètres, pas de la variance des données. La régularisation n’affecte pas la dispersion des données; elle réduit l’incertitude de nos estimés en les « tirant » vers zéro.
Spectre des valeurs singulières et rang effectif¶
En pratique, les données réelles ont souvent des dizaines ou des centaines de dimensions. Comment se comportent les valeurs singulières dans ce cas? Examinons un exemple avec des données de dimension plus élevée.
Source
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(123)
# Simuler des données avec une structure de rang bas + bruit
n_samples = 100
n_features = 30
# Vraie structure: combinaison de 5 facteurs latents
n_latent = 5
latent = np.random.randn(n_samples, n_latent)
loadings = np.random.randn(n_latent, n_features)
X_signal = latent @ loadings
# Ajouter du bruit
noise_level = 0.5
X = X_signal + noise_level * np.random.randn(n_samples, n_features)
# Centrer
X_centered = X - X.mean(axis=0)
# SVD
U, d, Vt = np.linalg.svd(X_centered, full_matrices=False)
# Variance expliquée
variance_explained = d**2 / np.sum(d**2)
cumulative_variance = np.cumsum(variance_explained)
# Créer la figure
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
# Panneau gauche: spectre des valeurs singulières
ax1 = axes[0]
ax1.semilogy(range(1, len(d)+1), d, 'o-', markersize=6, linewidth=1.5, color='tab:blue')
ax1.axhline(d[n_latent], color='tab:red', linestyle='--', linewidth=1.5,
label=f'Seuil (rang effectif = {n_latent})')
ax1.fill_between(range(1, n_latent+1), d[:n_latent], alpha=0.3, color='tab:green', label='Signal')
ax1.fill_between(range(n_latent+1, len(d)+1), d[n_latent:], alpha=0.3, color='tab:orange', label='Bruit')
ax1.set_xlabel('Indice $j$')
ax1.set_ylabel('Valeur singulière $d_j$ (échelle log)')
ax1.set_title('Spectre des valeurs singulières')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim(0.5, len(d)+0.5)
# Panneau droit: variance expliquée cumulative
ax2 = axes[1]
ax2.plot(range(1, len(d)+1), cumulative_variance * 100, 'o-', markersize=6,
linewidth=1.5, color='tab:blue')
ax2.axhline(95, color='gray', linestyle='--', linewidth=1, label='Seuil 95%')
ax2.axvline(n_latent, color='tab:red', linestyle='--', linewidth=1.5)
# Trouver k pour 95% de variance
k_95 = np.searchsorted(cumulative_variance, 0.95) + 1
ax2.scatter([k_95], [cumulative_variance[k_95-1]*100], s=100, color='tab:red', zorder=5)
ax2.annotate(f'{k_95} composantes\npour 95%', xy=(k_95, cumulative_variance[k_95-1]*100),
xytext=(k_95+5, cumulative_variance[k_95-1]*100-10), fontsize=10,
arrowprops=dict(arrowstyle='->', color='gray'))
ax2.set_xlabel('Nombre de composantes $k$')
ax2.set_ylabel('Variance expliquée cumulative (%)')
ax2.set_title('Variance expliquée')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0.5, len(d)+0.5)
ax2.set_ylim(0, 105)
plt.tight_layout()
Le panneau de gauche montre le spectre des valeurs singulières en échelle logarithmique. On observe un schéma typique:
Les premières valeurs singulières sont grandes: elles correspondent aux directions du signal, la vraie structure sous-jacente des données.
Après un certain point (ici, autour de ), les valeurs singulières chutent et forment un “plancher”: ce sont les directions du bruit.
Le rang effectif est le nombre de valeurs singulières significativement au-dessus du plancher de bruit. Dans cet exemple, nous avons simulé 5 facteurs latents, et le spectre révèle bien cette structure: les 5 premières valeurs singulières dominent.
Le panneau de droite montre la variance expliquée cumulative. C’est un outil pratique pour choisir combien de composantes retenir:
Critère du seuil: Retenir assez de composantes pour expliquer 95% (ou 99%) de la variance.
Critère du coude: Chercher le “coude” dans le spectre où les valeurs singulières cessent de décroître rapidement.
Critère du gap: Si le spectre présente un saut net (comme ici entre et ), c’est un bon point de coupure.
De la troncature à la réduction de dimension (ACP)¶
L’analyse ci-dessus suggère une idée: si les dernières directions ne contiennent que du bruit, pourquoi ne pas simplement les ignorer? Au lieu de rétrécir les coefficients comme Ridge, nous pourrions tronquer la représentation en ne gardant que les premières directions principales.
C’est exactement l’idée de l’analyse en composantes principales (ACP). Au lieu de travailler avec les caractéristiques originales, nous projetons les données sur les premiers vecteurs singuliers:
où contient les premiers vecteurs singuliers (les colonnes de correspondant aux plus grandes valeurs singulières).
Cette projection préserve au mieux la variance des données: les composantes principales capturent la direction où les données varient le plus. La reconstruction à partir de cette représentation compressée s’écrit:
L’erreur de reconstruction est minimale parmi toutes les projections linéaires sur un sous-espace de dimension .
Lien entre Ridge et ACP: Les deux approches traitent le même problème (les directions à faible valeur singulière sont bruitées) mais différemment:
| Approche | Traitement des directions bruitées | Type de régularisation |
|---|---|---|
| Ridge | Rétrécit (soft thresholding) | Continue: garde tout, pénalise |
| ACP | Élimine (hard thresholding) | Discrète: garde , ignore le reste |
Ridge est appropriée pour la régression supervisée, où même les petites directions peuvent contenir du signal utile pour prédire . L’ACP est appropriée pour la réduction de dimension non supervisée, où nous voulons une représentation compacte des données elles-mêmes.
Pourquoi aide¶
Ce terme diagonal a plusieurs effets bénéfiques:
Amélioration du conditionnement: La matrice peut être mal conditionnée (ses valeurs propres varient sur plusieurs ordres de grandeur) ou même singulière. L’ajout de augmente toutes les valeurs propres de , rendant la matrice inversible et mieux conditionnée.
Rétrécissement des coefficients (shrinkage): Comme nous l’avons vu dans la section SVD ci-dessus, la solution Ridge s’écrit:
Le facteur de rétrécissement est toujours inférieur à 1, ce qui “rétrécit” chaque composante vers zéro. L’effet est différencié selon les directions:
Pour une grande valeur singulière (fort signal), le facteur reste proche de 1 même pour des valeurs modérées de . La direction est peu affectée.
Pour une petite valeur singulière (faible signal), le facteur décroît rapidement avec . La direction est fortement pénalisée.
Pour visualiser ce rétrécissement et comprendre son effet, examinons un exemple concret. L’animation suivante montre simultanément trois perspectives sur la régularisation Ridge: les données et la droite ajustée, le paysage de perte avec la contrainte, et les facteurs de rétrécissement.
Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.patches import Circle
from IPython.display import Image
# Générer des données de régression simple
np.random.seed(42)
n = 30
# Une seule caractéristique pour visualisation claire
x = np.random.uniform(-2, 2, n)
# Relation linéaire avec bruit
theta_true = 1.5
y = theta_true * x + np.random.normal(0, 0.8, n)
# Ajouter une caractéristique corrélée (pour créer de la colinéarité)
x2 = 0.9 * x + 0.3 * np.random.randn(n)
# Matrice de design avec les deux caractéristiques
X = np.column_stack([x, x2])
# Solution MCO
theta_ols = np.linalg.lstsq(X, y, rcond=None)[0]
# SVD pour analyse
U, d_svd, Vt = np.linalg.svd(X, full_matrices=False)
V = Vt.T
# Fonction pour calculer la solution Ridge
def ridge_solution(X, y, lam):
n_features = X.shape[1]
return np.linalg.solve(X.T @ X + lam * np.eye(n_features), X.T @ y)
# Préparer la grille pour les contours RSS
theta1_range = np.linspace(-0.5, 3, 100)
theta2_range = np.linspace(-1.5, 2, 100)
T1, T2 = np.meshgrid(theta1_range, theta2_range)
# Calculer RSS pour chaque point de la grille
RSS = np.zeros_like(T1)
for i in range(T1.shape[0]):
for j in range(T1.shape[1]):
theta = np.array([T1[i, j], T2[i, j]])
residuals = y - X @ theta
RSS[i, j] = np.sum(residuals**2)
# Créer la figure avec trois panneaux
fig = plt.figure(figsize=(15, 5))
# === Panneau 1: Données et droite ajustée ===
ax1 = fig.add_subplot(1, 3, 1)
# Données
ax1.scatter(x, y, c='tab:blue', s=50, alpha=0.7, label='Données', zorder=3)
# Grille pour tracer les droites
x_grid = np.linspace(-2.5, 2.5, 100)
# Droite MCO (fixe) - on utilise seulement theta1 car x et x2 sont très corrélés
# La prédiction effective est environ (theta1 + 0.9*theta2) * x
slope_ols = theta_ols[0] + 0.9 * theta_ols[1] # Pente effective
y_ols = slope_ols * x_grid
ax1.plot(x_grid, y_ols, 'k-', linewidth=2, alpha=0.7, label='MCO')
# Droite Ridge (animée)
line_ridge, = ax1.plot([], [], '-', color='tab:orange', linewidth=2.5, label='Ridge')
# Ligne horizontale (prédiction = moyenne, lambda infini)
y_mean = np.mean(y)
ax1.axhline(y_mean, color='gray', linestyle=':', alpha=0.5, label=f'Moyenne ($\\lambda \\to \\infty$)')
ax1.set_xlabel('$x$')
ax1.set_ylabel('$y$')
ax1.set_title('Données et droite de régression')
ax1.legend(loc='upper left', fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-2.5, 2.5)
ax1.set_ylim(-4, 5)
# Texte pour les coefficients
coef_text = ax1.text(0.98, 0.02, '', transform=ax1.transAxes, fontsize=10,
ha='right', va='bottom',
bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
# === Panneau 2: Paysage de perte ===
ax2 = fig.add_subplot(1, 3, 2)
# Contours RSS (ellipses centrées sur OLS)
levels = np.percentile(RSS.flatten(), [5, 15, 30, 50, 70, 85, 95])
contours = ax2.contour(T1, T2, RSS, levels=levels, colors='gray', alpha=0.6)
ax2.clabel(contours, inline=True, fontsize=8, fmt='%.0f')
# Solution MCO (fixe)
ax2.plot(theta_ols[0], theta_ols[1], 'ko', markersize=12, label='MCO', zorder=5)
ax2.annotate('MCO', xy=(theta_ols[0], theta_ols[1]),
xytext=(theta_ols[0] + 0.2, theta_ols[1] + 0.2),
fontsize=11, ha='left')
# Origine = coefficients nuls (prédiction constante)
ax2.plot(0, 0, 'k+', markersize=15, markeredgewidth=2, zorder=4)
ax2.annotate('$\\boldsymbol{\\theta} = 0$\n(pente nulle)', xy=(0, 0),
xytext=(-0.4, -1.2), fontsize=9, ha='center', color='gray')
# Cercle de contrainte Ridge (animé)
circle_ridge = Circle((0, 0), radius=np.linalg.norm(theta_ols),
fill=False, edgecolor='tab:orange', linewidth=2.5,
linestyle='-', alpha=0.8, zorder=3)
ax2.add_patch(circle_ridge)
# Solution Ridge (animée)
point_ridge, = ax2.plot([], [], 'o', color='tab:orange', markersize=10,
label='Ridge', zorder=6)
# Chemin de régularisation
lambda_path = np.logspace(-3, 1.5, 50)
theta_path = np.array([ridge_solution(X, y, l) for l in lambda_path])
ax2.plot(theta_path[:, 0], theta_path[:, 1], 'tab:orange', linewidth=1.5,
alpha=0.4, linestyle='--', label='Chemin')
ax2.set_xlabel('$\\theta_1$')
ax2.set_ylabel('$\\theta_2$')
ax2.set_title('Paysage RSS et chemin de régularisation')
ax2.legend(loc='upper right', fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-0.5, 3)
ax2.set_ylim(-1.5, 2)
ax2.set_aspect('equal')
# Texte pour lambda
lambda_text = ax2.text(0.02, 0.98, '', transform=ax2.transAxes, fontsize=11,
va='top', ha='left',
bbox=dict(boxstyle='round', facecolor='white', alpha=0.9))
# === Panneau 3: Facteurs de rétrécissement ===
ax3 = fig.add_subplot(1, 3, 3)
lambda_range = np.linspace(0, 10, 200)
shrink1_curve = d_svd[0]**2 / (d_svd[0]**2 + lambda_range)
shrink2_curve = d_svd[1]**2 / (d_svd[1]**2 + lambda_range)
ax3.plot(lambda_range, shrink1_curve, 'b-', linewidth=2,
label=f'Direction forte ($d_1={d_svd[0]:.1f}$)')
ax3.plot(lambda_range, shrink2_curve, 'r-', linewidth=2,
label=f'Direction faible ($d_2={d_svd[1]:.2f}$)')
# Zone de surapprentissage et sous-apprentissage
ax3.axvspan(0, 0.5, alpha=0.1, color='red', label='Surapprentissage')
ax3.axvspan(5, 10, alpha=0.1, color='blue', label='Sous-apprentissage')
point_shrink1, = ax3.plot([], [], 'bo', markersize=10, zorder=3)
point_shrink2, = ax3.plot([], [], 'ro', markersize=10, zorder=3)
ax3.axhline(1.0, color='gray', linestyle='--', alpha=0.5)
ax3.set_xlabel('$\\lambda$')
ax3.set_ylabel('Facteur de rétrécissement')
ax3.set_title('Rétrécissement par direction SVD')
ax3.legend(loc='center right', fontsize=8)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 10)
ax3.set_ylim(0, 1.1)
shrink_text = ax3.text(0.02, 0.5, '', transform=ax3.transAxes, fontsize=10,
va='center', ha='left',
bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
plt.tight_layout()
# Fonction d'animation
def animate(frame):
if frame < 80:
lam = (frame / 80) * 10
else:
lam = 10
# Solution Ridge
theta_ridge = ridge_solution(X, y, lam)
# Panneau 1: Mettre à jour la droite
slope_ridge = theta_ridge[0] + 0.9 * theta_ridge[1]
y_ridge = slope_ridge * x_grid
line_ridge.set_data(x_grid, y_ridge)
coef_text.set_text(f'Pente MCO: {slope_ols:.2f}\nPente Ridge: {slope_ridge:.2f}')
# Panneau 2: Mettre à jour le cercle et le point
norm_ridge = np.linalg.norm(theta_ridge)
circle_ridge.set_radius(norm_ridge)
point_ridge.set_data([theta_ridge[0]], [theta_ridge[1]])
lambda_text.set_text(f'$\\lambda = {lam:.1f}$')
# Panneau 3: Mettre à jour les points de rétrécissement
shrink1 = d_svd[0]**2 / (d_svd[0]**2 + lam)
shrink2 = d_svd[1]**2 / (d_svd[1]**2 + lam)
point_shrink1.set_data([lam], [shrink1])
point_shrink2.set_data([lam], [shrink2])
shrink_text.set_text(f'Facteur dir. 1: {shrink1:.2f}\nFacteur dir. 2: {shrink2:.2f}')
return (line_ridge, point_ridge, circle_ridge, lambda_text,
point_shrink1, point_shrink2, coef_text, shrink_text)
# Créer l'animation
anim = FuncAnimation(fig, animate, frames=90, interval=80, blit=False, repeat=True)
anim.save('_static/ridge_geometry.gif', writer='pillow', fps=12, dpi=100)
plt.close()
# Afficher le GIF
Image(filename='_static/ridge_geometry.gif')
L’animation relie trois perspectives sur la régularisation Ridge lorsque augmente de 0 à 10:
Panneau de gauche (données et ajustement): Les points bleus sont les données d’entraînement. La droite noire est l’ajustement MCO (), la droite orange est l’ajustement Ridge. À mesure que augmente, la pente de la droite Ridge diminue, se rapprochant de la ligne horizontale (prédiction constante égale à la moyenne). C’est le rétrécissement vers zéro: Ridge “tire” les coefficients vers l’origine, ce qui réduit la pente.
Panneau central (paysage de perte): Chaque point de ce plan représente un choix de coefficients . Les contours gris montrent la fonction de coût RSS: plus on est proche du point noir (MCO), plus l’erreur sur les données d’entraînement est faible. L’ellipse est allongée car et sont corrélées (colinéarité). L’origine correspond à une pente nulle (prédiction constante). Le cercle orange montre la norme de la solution Ridge courante . La formulation pénalisée est équivalente à la formulation contrainte sous , où joue le rôle du multiplicateur de Lagrange: pour chaque , il existe un tel que les deux problèmes ont la même solution. À mesure que augmente, la solution se déplace le long du chemin de régularisation vers l’origine.
Panneau de droite (rétrécissement différencié): La direction “forte” (grande valeur singulière , où les données sont dispersées) est peu affectée par la régularisation. La direction “faible” (petite valeur singulière , direction de colinéarité) est rétrécit beaucoup plus rapidement. C’est le cœur de l’effet Ridge: pénaliser davantage les directions où le signal est faible et l’estimation instable.
L’intuition géométrique est la suivante: quand les données sont colinéaires, l’ellipse RSS est très allongée. De petites perturbations dans les données causent de grands déplacements de la solution MCO le long de l’axe allongé. La pénalité Ridge ajoute un terme qui « tire » la solution vers l’origine. Dans la formulation contrainte équivalente, cela correspond à chercher le minimum de RSS à l’intérieur d’une boule de rayon . Plus la boule est petite (plus est grand), plus la solution est proche de l’origine et donc plus stable.
Stabilité numérique: Quand est presque singulière, de petites perturbations dans les données causent de grandes variations dans . La régularisation réduit cette sensibilité.
Classes de modèles et expansion de caractéristiques¶
Les exemples de régularisation ci-dessus utilisaient des caractéristiques polynomiales: au lieu de prédire directement à partir de , nous avons construit des caractéristiques et appliqué un modèle linéaire dans cet espace étendu. Cette technique s’appelle l’expansion de caractéristiques et mérite d’être formalisée.
Trois familles de modèles¶
Situons les modèles linéaires dans une hiérarchie plus large. Nous distinguons trois familles de complexité croissante:
Modèles linéaires: . La sortie est une combinaison linéaire des entrées. Simple, interprétable, mais limité aux relations linéaires.
Modèles à expansion de caractéristiques: , où est une transformation non linéaire fixée à l’avance (par exemple, polynomiale). Le modèle reste linéaire dans les paramètres , ce qui facilite l’optimisation, mais peut capturer des relations non linéaires en . L’espace de redescription a souvent une dimension .
Réseaux de neurones: . Une composition de fonctions non linéaires, chacune avec ses propres paramètres. Contrairement aux modèles à expansion fixe, les réseaux de neurones apprennent la représentation en même temps que les paramètres .
Cette progression capture l’évolution historique du domaine: des modèles linéaires classiques aux méthodes à noyaux (expansion implicite), puis aux réseaux profonds qui apprennent leurs propres représentations. Nous verrons les réseaux de neurones en détail dans les chapitres suivants; concentrons-nous ici sur les deux premières familles.
Expansion de caractéristiques¶
Pour capturer des relations non linéaires tout en gardant un modèle linéaire dans les paramètres, nous transformons les entrées. En régression polynomiale, nous appliquons une fonction :
La prédiction devient . Le modèle est polynomial en mais linéaire en , ce qui permet d’utiliser les mêmes algorithmes d’optimisation (MCO, Ridge).
Le degré contrôle la capacité du modèle. Avec , nous avons une droite. Avec élevé, le polynôme peut osciller pour passer par tous les points d’entraînement. Avec , nous pouvons interpoler exactement les points: le risque empirique atteint zéro. Mais un polynôme qui passe exactement par les points d’entraînement n’a aucune raison de bien prédire les nouveaux points.
Illustrons ce phénomène avec les données de freinage. Nous ajustons des polynômes de degrés 1, 2, 5 et 15, et comparons leurs erreurs sur les ensembles d’entraînement et de test.
Source
import numpy as np
import matplotlib.pyplot as plt
import warnings
# Suppress polyfit warnings for high-degree polynomials (expected for this demo)
warnings.filterwarnings('ignore', message='Polyfit may be poorly conditioned')
# Données de freinage
speed = np.array([4, 4, 7, 7, 8, 9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14,
14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19,
20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], dtype=float)
dist = np.array([2, 10, 4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46,
26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46,
68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], dtype=float)
# Train/test split
np.random.seed(42)
indices = np.random.permutation(len(speed))
train_idx, test_idx = indices[:35], indices[35:]
speed_train, dist_train = speed[train_idx], dist[train_idx]
speed_test, dist_test = speed[test_idx], dist[test_idx]
degrees_to_plot = [1, 2, 5, 15]
degrees_eval = range(1, 16)
fig, axes = plt.subplots(2, 2, figsize=(10, 8))
# Pre-compute all errors for the summary plot later
all_train_errors = []
all_test_errors = []
for deg in degrees_eval:
coeffs = np.polyfit(speed_train, dist_train, deg)
all_train_errors.append(np.mean((dist_train - np.polyval(coeffs, speed_train))**2))
all_test_errors.append(np.mean((dist_test - np.polyval(coeffs, speed_test))**2))
for ax, deg in zip(axes.flat, degrees_to_plot):
# Fit polynomial
coeffs = np.polyfit(speed_train, dist_train, deg)
# Predictions
pred_train = np.polyval(coeffs, speed_train)
pred_test = np.polyval(coeffs, speed_test)
# MSE
mse_train = np.mean((dist_train - pred_train)**2)
mse_test = np.mean((dist_test - pred_test)**2)
# Plot
ax.scatter(speed_train, dist_train, alpha=0.6, s=30, label='Entraînement')
ax.scatter(speed_test, dist_test, alpha=0.6, s=30, marker='s', label='Test')
speed_grid = np.linspace(3, 26, 200)
pred_grid = np.polyval(coeffs, speed_grid)
# Clip extreme predictions for visualization
pred_grid = np.clip(pred_grid, -50, 200)
ax.plot(speed_grid, pred_grid, 'k-', alpha=0.7)
ax.set_xlim(3, 26)
ax.set_ylim(-20, 150)
ax.set_xlabel('Vitesse (mph)')
ax.set_ylabel('Distance (ft)')
ax.set_title(f'Degré {deg}: Entr. EQM={mse_train:.1f}, Test EQM={mse_test:.1f}')
if deg == 1:
ax.legend()
plt.tight_layout()
Le polynôme de degré 1 (droite) ne capture pas la courbure des données: c’est du sous-apprentissage. Le polynôme de degré 2 capture bien la relation quadratique. Le polynôme de degré 5 commence à osciller. Le polynôme de degré 15 passe près de tous les points d’entraînement, mais ses oscillations produisent des prédictions absurdes entre les points: c’est du surapprentissage.
Source
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(degrees_eval, all_train_errors, 'o-', linewidth=2, label='Erreur entraînement')
ax.plot(degrees_eval, all_test_errors, 's-', linewidth=2, label='Erreur test')
# Utiliser une échelle logarithmique car l'erreur de test explose
ax.set_yscale('log')
ax.set_xlabel('Degré du polynôme (complexité)')
ax.set_ylabel('EQM (échelle log)')
ax.set_xticks(range(1, 16, 2))
ax.grid(True, which="both", ls="-", alpha=0.2)
ax.legend()
ax.set_title('Compromis biais-variance')
plt.tight_layout()
L’erreur d’entraînement diminue avec le degré du polynôme. L’erreur de test diminue d’abord (quand le modèle gagne en expressivité), puis augmente (quand le modèle commence à mémoriser le bruit). Le meilleur modèle se trouve à l’intersection de ces deux tendances. La régularisation Ridge, vue précédemment, est une alternative au choix du degré: elle permet d’utiliser un modèle de haute capacité tout en contrôlant le surapprentissage.
Décomposition biais-variance¶
Ce compromis peut être formalisé mathématiquement. Supposons que les données suivent le modèle , où est la vraie fonction (le prédicteur de Bayes optimal pour la perte quadratique), et est un bruit de moyenne nulle et de variance .
Clarification sur les variables aléatoires: Dans cette analyse, il y a trois sources d’aléa qu’il faut bien distinguer:
est une valeur déterministe (fixe) pour chaque . C’est la vraie fonction sous-jacente, qui ne dépend d’aucun tirage aléatoire.
est une variable aléatoire représentant le bruit d’observation. C’est l’aléa intrinsèque aux données.
est notre estimateur, une fonction apprise à partir de l’échantillon d’entraînement . Comme est tiré aléatoirement, différents tirages de produisent différentes fonctions . Pour un point test fixé, la prédiction est donc une variable aléatoire (un scalaire qui varie selon ), même si elle-même est une fonction.
Autrement dit, l’espérance moyenne sur tous les échantillons d’entraînement possibles: si nous pouvions répéter l’expérience d’apprentissage un grand nombre de fois avec différents , quelle serait la prédiction moyenne pour ce ?
L’erreur quadratique moyenne, moyennée sur les échantillons d’entraînement possibles et le bruit, se décompose comme suit:
En développant le carré et utilisant ainsi que l’indépendance entre et :
Pour le premier terme, ajoutons et retranchons :
En développant et utilisant :
Nous obtenons la décomposition biais-variance:
Chaque terme a une interprétation précise:
Biais²: L’écart entre la prédiction moyenne de notre estimateur et la vraie fonction . Un modèle trop simple (classe trop restrictive) aura un biais élevé car il ne peut pas approcher .
Variance: La sensibilité de notre estimateur à l’échantillon d’entraînement particulier. Un modèle trop complexe aura une variance élevée car de petites variations dans les données causent de grandes variations dans les prédictions.
: Le bruit irréductible, inhérent aux données. Aucun estimateur ne peut faire mieux que cette erreur.
Cette décomposition explique pourquoi l’erreur de test a une forme en U: à faible complexité, le biais domine; à haute complexité, la variance domine. Le minimum se trouve au point où la somme des deux est minimale.
Intuition géométrique: pourquoi la dimension supérieure aide¶
L’expansion de caractéristiques semble être un simple changement de variables, mais elle cache une idée géométrique profonde. Pour comprendre pourquoi projeter les données dans un espace de dimension supérieure permet de capturer des relations non linéaires, examinons d’abord le cas de la régression, puis celui de la classification.
Le plan caché derrière la parabole¶
Considérons une régression quadratique: . Cette fonction est non linéaire en (c’est une parabole), mais linéaire dans les paramètres . Que signifie cette distinction géométriquement?
Introduisons l’espace des caractéristiques . Chaque valeur de correspond à un point dans . Ces points ne sont pas dispersés arbitrairement: ils vivent sur une courbe particulière, la courbe des moments (moment curve), qui ressemble à une rampe tordue.
Source
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Créer la courbe des moments: (x, x², x³) pour visualisation
# On utilise (1, x, x²) mais on projette sur (x, x², y) pour la visualisation
fig = plt.figure(figsize=(12, 5))
# Gauche: les fonctions de base
ax1 = fig.add_subplot(121)
x = np.linspace(-2, 2, 100)
ax1.plot(x, np.ones_like(x), 'b-', linewidth=2, label=r'$\phi_0(x) = 1$')
ax1.plot(x, x, 'orange', linewidth=2, label=r'$\phi_1(x) = x$')
ax1.plot(x, x**2, 'g-', linewidth=2, label=r'$\phi_2(x) = x^2$')
ax1.axhline(0, color='gray', linewidth=0.5)
ax1.axvline(0, color='gray', linewidth=0.5)
ax1.set_xlabel('$x$')
ax1.set_ylabel(r'$\phi_j(x)$')
ax1.set_title('Les fonctions de base')
ax1.legend()
ax1.set_ylim(-2.5, 4.5)
ax1.grid(True, alpha=0.3)
# Droite: combinaisons linéaires
ax2 = fig.add_subplot(122)
x = np.linspace(-2, 2, 100)
# Différentes combinaisons de coefficients
combinations = [
((1, 0, 0), 'Constante: $1$'),
((0, 1, 0), 'Linéaire: $x$'),
((0, 0, 1), 'Quadratique: $x^2$'),
((1, -0.5, 0.5), 'Combinaison: $1 - 0.5x + 0.5x^2$'),
]
colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red']
for (theta0, theta1, theta2), label in combinations:
y = theta0 + theta1 * x + theta2 * x**2
ax2.plot(x, y, linewidth=2, label=label)
ax2.axhline(0, color='gray', linewidth=0.5)
ax2.axvline(0, color='gray', linewidth=0.5)
ax2.set_xlabel('$x$')
ax2.set_ylabel('$f(x)$')
ax2.set_title('Combinaisons linéaires des fonctions de base')
ax2.legend(loc='upper center')
ax2.set_ylim(-2.5, 4.5)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
Les fonctions de base sont les “ingrédients” du modèle. La régression polynomiale cherche les coefficients qui mélangent ces ingrédients de façon optimale. Chaque combinaison produit une courbe différente, mais toutes sont des paraboles (ou des cas dégénérés: droites, constantes).
Voici l’insight géométrique clé: dans l’espace , le modèle définit un plan. La parabole que nous voyons dans le graphique est simplement la projection de ce plan sur notre espace de visualisation.
Source
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from IPython.display import Image
# Générer des données avec relation quadratique
np.random.seed(42)
n_points = 40
x_data = np.random.uniform(-1.8, 1.8, n_points)
y_true = 0.5 + 0.3 * x_data + 0.8 * x_data**2
y_data = y_true + np.random.normal(0, 0.3, n_points)
# Ajuster le modèle quadratique
X_design = np.column_stack([np.ones(n_points), x_data, x_data**2])
theta = np.linalg.lstsq(X_design, y_data, rcond=None)[0]
# Créer l'animation
fig = plt.figure(figsize=(9, 7))
ax = fig.add_subplot(111, projection='3d')
# Grille pour le plan de régression
x_grid = np.linspace(-2, 2, 20)
x2_grid = np.linspace(0, 4, 20)
X_plane, X2_plane = np.meshgrid(x_grid, x2_grid)
Y_plane = theta[0] + theta[1] * X_plane + theta[2] * X2_plane
# Surface parabolique z = x²
x_surf = np.linspace(-2, 2, 30)
X_surf, Y_surf_temp = np.meshgrid(x_surf, np.linspace(-1, 5, 30))
Z_surf = X_surf**2
def init():
ax.clear()
return []
def animate(frame):
ax.clear()
# Animation en trois phases
if frame < 15:
# Phase 1: Vue 2D (de côté, cachant x²)
elev = 0
azim = 0
show_plane = False
show_surface = False
title = 'Vue 2D: régression quadratique'
elif frame < 35:
# Phase 2: Rotation révélant la 3ème dimension
progress = (frame - 15) / 20
elev = progress * 25
azim = progress * 45
show_plane = False
show_surface = True
title = 'Rotation: découverte de la dimension $x^2$...'
else:
# Phase 3: Vue 3D avec le plan
elev = 25
azim = 45 + (frame - 35) * 2
show_plane = True
show_surface = True
title = r'Le plan $y = \theta_0 + \theta_1 x + \theta_2 x^2$ dans lespace 3D'
ax.view_init(elev=elev, azim=azim)
# Surface parabolique (rampe z = x²)
if show_surface:
ax.plot_surface(X_surf, Z_surf, Y_surf_temp, alpha=0.1, color='gray')
# Points de données dans l'espace 3D: (x, x², y)
ax.scatter(x_data, x_data**2, y_data, c='tab:blue', s=50, alpha=0.8,
label='Données', depthshade=True)
# Plan de régression
if show_plane:
ax.plot_surface(X_plane, X2_plane, Y_plane, alpha=0.4, color='tab:orange',
label='Plan de régression')
# Courbe de régression sur la surface z = x²
x_curve = np.linspace(-1.8, 1.8, 100)
y_curve = theta[0] + theta[1] * x_curve + theta[2] * x_curve**2
ax.plot(x_curve, x_curve**2, y_curve, 'r-', linewidth=3,
label='Courbe ajustée')
ax.set_xlabel('$x$')
ax.set_ylabel('$x^2$')
ax.set_zlabel('$y$')
ax.set_xlim(-2, 2)
ax.set_ylim(0, 4)
ax.set_zlim(-1, 5)
ax.set_title(title, fontsize=11)
return []
anim = FuncAnimation(fig, animate, init_func=init, frames=70, interval=100, blit=True)
anim.save('_static/regression_plane_3d.gif', writer='pillow', fps=10, dpi=100)
plt.close()
# Afficher le GIF
Image(filename='_static/regression_plane_3d.gif')
L’animation révèle la structure cachée de la régression polynomiale:
Vue 2D initiale: On voit les données et la parabole ajustée, comme dans un graphique classique.
Rotation: En faisant pivoter la vue, on découvre que chaque point vit en réalité dans un espace 3D aux coordonnées . La surface grise représente la “rampe” sur laquelle tous les points sont contraints de vivre.
Le plan de régression: Le modèle est un plan (en orange) dans cet espace 3D. Ce plan est choisi pour minimiser les distances verticales aux points.
La courbe ajustée: La parabole rouge est l’intersection du plan avec la rampe . C’est ce que nous voyons quand nous projetons le tout sur le plan .
Cette perspective unifie le “linéaire dans les paramètres” et le “non linéaire en ”:
Linéaire dans les paramètres: Le modèle est un plan (objet linéaire) dans l’espace des caractéristiques.
Non linéaire en : La contrainte force les données à vivre sur une surface courbe, et l’intersection du plan avec cette surface produit une courbe.
De la régression à la classification¶
La même intuition s’applique en classification, avec une différence: au lieu de chercher un plan qui ajuste les données, on cherche un hyperplan qui les sépare. Voyons comment l’expansion de caractéristiques transforme des données non linéairement séparables en données linéairement séparables.
De 1D à 2D: séparer l’inséparable¶
Considérons des points sur une droite, répartis en deux classes: les points bleus au centre, les points orange aux extrémités. Aucun seuil unique ne peut séparer ces deux classes.
Source
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# Classe bleue: points au centre
x_blue = np.random.uniform(-0.5, 0.5, 15)
# Classe orange: points aux extrémités
x_orange = np.concatenate([np.random.uniform(-1.5, -0.8, 8),
np.random.uniform(0.8, 1.5, 8)])
fig, ax = plt.subplots(figsize=(10, 2))
ax.scatter(x_blue, np.zeros_like(x_blue), c='tab:blue', s=80, label='Classe A', zorder=3)
ax.scatter(x_orange, np.zeros_like(x_orange), c='tab:orange', s=80, label='Classe B', zorder=3)
ax.axhline(0, color='gray', linewidth=0.5, zorder=1)
ax.set_xlim(-2, 2)
ax.set_ylim(-0.5, 0.5)
ax.set_xlabel('$x$')
ax.set_yticks([])
ax.legend(loc='upper right')
ax.set_title('Données 1D: aucun seuil ne sépare les deux classes')
plt.tight_layout()
Appliquons maintenant l’expansion . Chaque point est projeté sur une parabole dans l’espace 2D. Les points proches de zéro (classe bleue) restent bas, tandis que les points éloignés de zéro (classe orange) montent.
Source
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# Gauche: espace original avec tentative de séparation
ax = axes[0]
ax.scatter(x_blue, np.zeros_like(x_blue), c='tab:blue', s=80, zorder=3)
ax.scatter(x_orange, np.zeros_like(x_orange), c='tab:orange', s=80, zorder=3)
ax.axhline(0, color='gray', linewidth=0.5, zorder=1)
ax.axvline(0.3, color='red', linestyle='--', linewidth=2, label='Seuil?')
ax.set_xlim(-2, 2)
ax.set_ylim(-0.5, 0.5)
ax.set_xlabel('$x$')
ax.set_yticks([])
ax.set_title('Espace original: pas de séparation linéaire')
ax.legend()
# Droite: espace transformé
ax = axes[1]
ax.scatter(x_blue, x_blue**2, c='tab:blue', s=80, zorder=3, label='Classe A')
ax.scatter(x_orange, x_orange**2, c='tab:orange', s=80, zorder=3, label='Classe B')
# Ligne de séparation dans l'espace transformé
x_line = np.linspace(-2, 2, 100)
threshold = 0.6
ax.axhline(threshold, color='green', linestyle='-', linewidth=2, label='Frontière linéaire')
ax.fill_between(x_line, 0, threshold, alpha=0.1, color='blue')
ax.fill_between(x_line, threshold, 2.5, alpha=0.1, color='orange')
# Parabole de référence
x_curve = np.linspace(-1.6, 1.6, 100)
ax.plot(x_curve, x_curve**2, 'k--', alpha=0.3, linewidth=1)
ax.set_xlim(-2, 2)
ax.set_ylim(-0.1, 2.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$x^2$')
ax.set_title(r'Espace transformé $\phi(x) = (x, x^2)$: séparation linéaire!')
ax.legend()
plt.tight_layout()
Dans l’espace transformé, une simple droite horizontale sépare les deux classes. Cette droite correspond, dans l’espace original, à deux seuils: , soit . L’expansion de caractéristiques a transformé une frontière de décision non linéaire (un intervalle) en une frontière linéaire (une droite).
De 2D à 3D: soulever pour séparer¶
Passons à un exemple plus visuel. Considérons deux classes disposées en cercles concentriques: la classe bleue forme un disque central, la classe orange forme un anneau extérieur. Aucune droite ne peut séparer ces deux régions.
Source
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(123)
# Classe bleue: disque central
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)
# Classe orange: anneau extérieur
n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)
fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(x_blue, y_blue, c='tab:blue', s=40, alpha=0.7, label='Classe A')
ax.scatter(x_orange, y_orange, c='tab:orange', s=40, alpha=0.7, label='Classe B')
# Montrer qu'une droite ne peut pas séparer
theta_line = np.pi/4
x_line = np.linspace(-2, 2, 100)
y_line = np.tan(theta_line) * x_line
ax.plot(x_line, y_line, 'r--', linewidth=2, alpha=0.7, label='Droite?')
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_aspect('equal')
ax.legend()
ax.set_title('Cercles concentriques: pas de séparation linéaire en 2D')
plt.tight_layout()
Appliquons l’expansion . La troisième coordonnée est le carré de la distance à l’origine: . Les points proches du centre (petit ) sont “soulevés” moins haut que les points éloignés (grand ).
L’animation suivante montre cette transformation. En faisant pivoter la vue, on voit que les données, une fois projetées dans l’espace 3D, deviennent séparables par un plan horizontal.
Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import Image
np.random.seed(123)
# Classe bleue: disque central
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)
z_blue = x_blue**2 + y_blue**2
# Classe orange: anneau extérieur
n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)
z_orange = x_orange**2 + y_orange**2
# Créer la figure
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
# Plan de séparation
z_sep = 0.75
xx, yy = np.meshgrid(np.linspace(-1.8, 1.8, 10), np.linspace(-1.8, 1.8, 10))
zz = np.ones_like(xx) * z_sep
def init():
ax.clear()
return []
def animate(frame):
ax.clear()
# Élévation: commence à 90° (vue de dessus), descend à 25°
if frame < 20:
elev = 90 - frame * 3.25 # 90 -> 25
azim = 45
else:
elev = 25
azim = 45 + (frame - 20) * 4 # rotation horizontale
ax.view_init(elev=elev, azim=azim)
# Points
ax.scatter(x_blue, y_blue, z_blue, c='tab:blue', s=30, alpha=0.8, label='Classe A')
ax.scatter(x_orange, y_orange, z_orange, c='tab:orange', s=30, alpha=0.8, label='Classe B')
# Plan de séparation (apparaît après la descente)
if frame >= 15:
alpha_plane = min(0.3, (frame - 15) * 0.06)
ax.plot_surface(xx, yy, zz, alpha=alpha_plane, color='green')
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$x_1^2 + x_2^2$')
ax.set_xlim(-1.8, 1.8)
ax.set_ylim(-1.8, 1.8)
ax.set_zlim(0, 2.5)
# Titre dynamique
if frame < 10:
ax.set_title('Vue de dessus: cercles concentriques', fontsize=11)
elif frame < 20:
ax.set_title('Rotation: on découvre la structure 3D...', fontsize=11)
else:
ax.set_title(r'Espace $\phi(x_1,x_2) = (x_1, x_2, x_1^2+x_2^2)$: un plan sépare!', fontsize=11)
return []
anim = FuncAnimation(fig, animate, init_func=init, frames=60, interval=100, blit=True)
anim.save('_static/feature_expansion_3d.gif', writer='pillow', fps=10, dpi=100)
plt.close()
# Afficher le GIF
Image(filename='_static/feature_expansion_3d.gif')
L’animation montre comment l’expansion de caractéristiques “déplie” la géométrie des données. Vue de dessus, la structure est celle des cercles concentriques originaux. Mais la troisième dimension sépare verticalement les deux classes: le disque central reste bas, l’anneau extérieur monte. Un plan horizontal (en vert) suffit alors à séparer les classes.
Ce plan horizontal correspond, dans l’espace original 2D, à un cercle de rayon . La frontière de décision linéaire en 3D devient une frontière circulaire en 2D.
Source
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(123)
# Recréer les données (même seed que l'animation)
n_blue = 50
r_blue = np.random.uniform(0, 0.7, n_blue)
theta_blue = np.random.uniform(0, 2*np.pi, n_blue)
x_blue = r_blue * np.cos(theta_blue)
y_blue = r_blue * np.sin(theta_blue)
n_orange = 70
r_orange = np.random.uniform(1.0, 1.5, n_orange)
theta_orange = np.random.uniform(0, 2*np.pi, n_orange)
x_orange = r_orange * np.cos(theta_orange)
y_orange = r_orange * np.sin(theta_orange)
fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(x_blue, y_blue, c='tab:blue', s=40, alpha=0.7, label='Classe A')
ax.scatter(x_orange, y_orange, c='tab:orange', s=40, alpha=0.7, label='Classe B')
# Cercle de décision (projection du plan z = 0.75)
theta_circle = np.linspace(0, 2*np.pi, 100)
r_decision = np.sqrt(0.75)
ax.plot(r_decision * np.cos(theta_circle), r_decision * np.sin(theta_circle),
'g-', linewidth=2.5, label=f'Frontière: $r = {r_decision:.2f}$')
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_aspect('equal')
ax.legend()
ax.set_title('Frontière de décision projetée en 2D')
plt.tight_layout()
Le principe unificateur¶
Ces exemples de régression et de classification illustrent le même principe géométrique:
| Régression | Classification | |
|---|---|---|
| Objectif | Ajuster les données | Séparer les classes |
| Dans l’espace original | Courbe (parabole, etc.) | Frontière courbe (cercle, etc.) |
| Dans l’espace des caractéristiques | Hyperplan d’ajustement | Hyperplan séparateur |
| La courbe/frontière est... | L’intersection du plan avec la surface | La projection de l’hyperplan |
L’expansion de caractéristiques transforme un problème non linéaire en un problème linéaire dans un espace de dimension supérieure. Les modèles linéaires, simples à optimiser et à analyser, deviennent alors suffisants pour capturer des structures complexes.
En augmentant la dimension de l’espace de représentation, nous augmentons la capacité du modèle. Pour la classification, un résultat classique de géométrie affirme que points en position générale sont presque sûrement séparables par un hyperplan si la dimension de l’espace est au moins . Pour la régression, un polynôme de degré peut interpoler exactement points.
Mais cette flexibilité a un coût: plus l’espace est grand, plus le modèle risque de mémoriser les particularités des données d’entraînement plutôt que d’apprendre la structure sous-jacente. C’est le compromis biais-variance, et c’est pourquoi la régularisation (vue précédemment) est si importante pour les modèles à haute capacité.
Évaluation et choix de modèle¶
En pratique, nous estimons le risque par le risque empirique sur un ensemble de test disjoint de l’ensemble d’entraînement. Un troisième ensemble, l’ensemble de validation, sert à choisir parmi plusieurs modèles ou à régler des hyperparamètres. L’ensemble de test doit rester intact jusqu’à l’évaluation finale, pour fournir une estimation non biaisée.
Cette séparation est importante. Si nous utilisons l’ensemble de test pour faire des choix (quel modèle garder, quelle valeur d’hyperparamètre utiliser), l’estimation de performance sur ce même ensemble devient optimiste.
Hyperparamètres et validation¶
De nombreux modèles ont des hyperparamètres: des choix qui doivent être faits avant l’entraînement et qui ne sont pas appris à partir des données. Le degré d’un polynôme, le nombre de voisins dans les plus proches voisins, ou le coefficient de régularisation sont des exemples d’hyperparamètres.
Un hyperparamètre mal choisi peut mener au surapprentissage (modèle trop complexe) ou au sous-apprentissage (modèle trop simple). La méthode standard pour choisir un hyperparamètre est la validation: on réserve une partie des données (typiquement 20%) comme ensemble de validation, on entraîne le modèle pour plusieurs valeurs de l’hyperparamètre, et on retient celle qui minimise l’erreur sur l’ensemble de validation.
Plus formellement, pour un hyperparamètre , définissons le risque de validation:
où est le modèle entraîné avec l’hyperparamètre . La recherche par grille consiste à évaluer ce risque pour un ensemble de valeurs candidates et à retenir:
Une fois choisi, on peut ré-entraîner le modèle sur l’ensemble des données (entraînement + validation) pour obtenir le modèle final.
Validation croisée¶
Quand les données sont peu nombreuses, réserver 20% pour la validation peut être coûteux. La validation croisée offre une alternative.
L’idée est de partitionner les données en blocs. Pour chaque bloc , on entraîne le modèle sur les autres blocs et on évalue sur le bloc . Le risque de validation croisée est la moyenne des évaluations:
où désigne toutes les données sauf le bloc .
Source
import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 3))
K = 5
for fold in range(K):
for k in range(K):
if k == fold:
ax.barh(fold, 1, left=k, color='C1', edgecolor='black', linewidth=1, label='Validation' if fold == 0 else '')
else:
ax.barh(fold, 1, left=k, color='C0', edgecolor='black', linewidth=1, label='Entraînement' if fold == 0 and k == 0 else '')
ax.set_yticks(range(K))
ax.set_yticklabels([f'Itération {k+1}' for k in range(K)])
ax.set_xticks(np.arange(K) + 0.5)
ax.set_xticklabels([f'Bloc {k+1}' for k in range(K)])
ax.set_xlabel('Blocs de données')
ax.legend(loc='upper right')
ax.set_title(f'Validation croisée à {K} blocs')
ax.set_xlim(0, K)
plt.tight_layout()
Le cas particulier (un bloc par exemple) est appelé validation croisée leave-one-out. Elle utilise au maximum les données disponibles, mais est coûteuse en calcul. En pratique, ou offrent un bon compromis.
Biais inductifs¶
Il n’existe pas de modèle universel qui fonctionne optimalement pour tous les problèmes. Ce résultat, connu sous le nom de théorème du no free lunch, affirme qu’un algorithme d’apprentissage qui performe bien sur une classe de problèmes performe nécessairement moins bien sur d’autres.
Tout modèle encode des biais inductifs: des hypothèses implicites ou explicites sur la structure du problème. La régression linéaire suppose que la relation entre entrées et sorties est linéaire. Les k plus proches voisins supposent que les points proches dans l’espace des entrées ont des sorties similaires. Les modèles plus complexes, comme les réseaux de neurones, encodent d’autres hypothèses sur la structure des données.
Ces hypothèses sont nécessaires pour que l’apprentissage soit possible. Sans elles, nous n’aurions aucune raison de croire que la performance sur l’échantillon d’entraînement prédit la performance sur de nouvelles données. Le choix du modèle et de ses hypothèses est une décision que l’algorithme ne peut pas prendre seul; elle requiert une connaissance du domaine.
Fonctions de perte de substitution¶
La perte 0-1 pose un problème pratique. Les méthodes d’optimisation itératives, comme la descente de gradient, requièrent que la fonction objectif soit différentiable. Or la perte 0-1 est constante par morceaux: sa dérivée est nulle presque partout et indéfinie aux points de discontinuité.
Nous contournons ce problème en utilisant des fonctions de perte de substitution: des approximations convexes et différentiables de la perte originale.
Pour la classification binaire, plutôt que de prédire directement une classe, les modèles produisent souvent un score (un nombre réel). La prédiction de classe se fait ensuite en prenant le signe de ce score: si , on prédit la classe +1; si , on prédit la classe -1. La valeur absolue de mesure la confiance: plus est grand, plus le modèle est confiant dans sa prédiction.
Pour la classification binaire avec , la perte logistique est:
où est le score produit par le modèle. Cette fonction est convexe et différentiable partout. Lorsque et ont le même signe (prédiction correcte avec confiance), la perte est faible. Lorsqu’ils ont des signes opposés (erreur), la perte croît linéairement avec l’amplitude de l’erreur.
La perte à charnière (hinge loss) est utilisée dans les machines à vecteurs de support:
Cette fonction est convexe mais non différentiable au point . Elle est nulle lorsque la prédiction est correcte avec une marge suffisante (), et croît linéairement sinon.
Ces deux fonctions majorent la perte 0-1: pour tout et , nous avons et . Minimiser ces substituts garantit donc un certain contrôle sur la perte originale.
Source
import numpy as np
import matplotlib.pyplot as plt
# Margin: y * s (positive = correct prediction, negative = error)
margin = np.linspace(-3, 3, 500)
# 0-1 loss: 1 if margin < 0, else 0
loss_01 = (margin < 0).astype(float)
# Logistic loss: log(1 + exp(-margin))
loss_log = np.log(1 + np.exp(-margin))
# Hinge loss: max(0, 1 - margin)
loss_hinge = np.maximum(0, 1 - margin)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(margin, loss_01, 'k-', linewidth=2, label='Perte 0-1')
ax.plot(margin, loss_log, 'C0-', linewidth=2, label='Perte logistique')
ax.plot(margin, loss_hinge, 'C1-', linewidth=2, label='Perte à charnière')
ax.axvline(0, color='gray', linestyle=':', alpha=0.5)
ax.axhline(1, color='gray', linestyle=':', alpha=0.3)
ax.set_xlabel(r'Marge $y \cdot s$')
ax.set_ylabel('Perte')
ax.set_xlim(-3, 3)
ax.set_ylim(-0.1, 4)
ax.legend()
ax.set_title('Fonctions de perte de substitution comme bornes supérieures convexes')
# Annotate regions
ax.text(-1.5, 3.5, 'Erreur\n(prédiction incorrecte)', ha='center', fontsize=9, color='gray')
ax.text(1.5, 0.3, 'Correct\n(prédiction juste)', ha='center', fontsize=9, color='gray')
plt.tight_layout()
La figure montre les trois fonctions de perte en fonction de la marge . Une marge positive indique une prédiction correcte (le signe de correspond à ), une marge négative indique une erreur. La perte 0-1 est discontinue au point . Les pertes logistique et à charnière sont continues et convexes, ce qui permet d’utiliser des méthodes d’optimisation par gradient. Elles majorent partout la perte 0-1.
Le cadre probabiliste¶
Jusqu’ici, nous avons choisi des fonctions de perte de manière ad hoc: la perte quadratique semble raisonnable pour la régression, la perte logistique pour la classification. Mais d’où viennent ces choix? Existe-t-il un principe unificateur?
Le cadre probabiliste offre une réponse: plutôt que de choisir une perte arbitraire, nous modélisons explicitement comment les données ont été générées. Cette section présente d’abord le cadre général de l’inférence bayésienne, puis développe deux approches concrètes: le maximum de vraisemblance (EMV) et le maximum a posteriori (MAP).
Le cadre bayésien¶
La statistique bayésienne propose un cadre général pour l’estimation de paramètres. Au lieu d’estimer un point unique, elle caractérise notre incertitude sur les paramètres par une distribution de probabilité.
Le théorème de Bayes nous dit comment mettre à jour nos croyances sur les paramètres après avoir observé des données :
Chaque terme a un nom et un rôle précis:
est la distribution a posteriori: notre croyance sur après avoir vu les données
est la distribution a priori: notre croyance sur avant d’observer les données
est la vraisemblance: la probabilité des données pour un choix de paramètres
est la vraisemblance marginale: une constante de normalisation
L’a priori encode notre connaissance préalable. Pour une pièce de monnaie, nous pourrions croire que est probablement proche de 0.5. L’a posteriori combine cette croyance avec l’évidence des données.
Le prédicteur de Bayes optimal¶
Commençons par un idéal théorique. Si nous connaissions la vraie distribution conjointe , quelle fonction minimiserait le risque?
Puisque est toujours positif, minimiser cette intégrale revient à minimiser, pour chaque , l’espérance conditionnelle de la perte. Le prédicteur de Bayes optimal est donc:
La réponse dépend de la fonction de perte choisie.
Pour la perte quadratique , développons l’espérance:
C’est une fonction quadratique en . En dérivant et en posant la dérivée égale à zéro:
Le prédicteur optimal est la moyenne conditionnelle.
Pour la perte 0-1 en classification , l’espérance est:
Minimiser cette quantité revient à maximiser , donc à choisir la classe la plus probable:
Le prédicteur optimal est le mode conditionnel. Chaque perte définit son propre prédicteur optimal.
Ce prédicteur est un repère théorique: aucun algorithme ne peut faire mieux, car il suppose l’accès à la vraie distribution. La différence entre le risque d’un prédicteur appris et ce risque de Bayes mesure ce que nous perdons en ne connaissant pas la vraie distribution.
Prédiction bayésienne et distribution prédictive a posteriori¶
En pratique, nous ne connaissons pas . Nous avons un modèle paramétrique et une distribution a posteriori sur les paramètres. L’approche bayésienne complète consiste à moyenner les prédictions sur tous les paramètres possibles, pondérés par leur probabilité a posteriori:
Cette distribution prédictive a posteriori intègre l’incertitude sur les paramètres. Elle ne s’engage pas sur une valeur unique de , mais considère toutes les valeurs plausibles.
Le problème: cette intégrale est rarement calculable analytiquement. Elle nécessite de sommer sur un espace de paramètres de grande dimension, ce qui est coûteux ou impossible en pratique. C’est pourquoi nous recourons souvent à des estimateurs ponctuels: plutôt que d’intégrer sur tous les , nous en choisissons un seul, comme l’EMV ou le MAP.
Utilité du modèle probabiliste¶
Si nous finissons souvent par utiliser un estimateur ponctuel, pourquoi adopter le cadre probabiliste? Plusieurs raisons:
Justifier la fonction de perte: La perte quadratique découle naturellement de l’hypothèse de bruit gaussien. La perte logarithmique vient du principe de maximum de vraisemblance. Le cadre probabiliste explique pourquoi ces choix sont raisonnables.
Quantifier l’incertitude: Au-delà de la prédiction ponctuelle , nous pouvons donner un intervalle de prédiction. Sous un modèle gaussien, a environ 95% de chances de tomber dans .
Comparer des modèles: La vraisemblance marginale permet de comparer des modèles de complexités différentes, pénalisant automatiquement les modèles trop complexes.
Ouvrir la porte à l’inférence complète: Quand les ressources le permettent (méthodes de Monte Carlo, inférence variationnelle), nous pouvons approximer la distribution prédictive complète plutôt que de nous limiter à un point.
Maximum de vraisemblance¶
Le maximum de vraisemblance est la première approche concrète dans le cadre probabiliste: nous cherchons les paramètres qui rendent nos observations les plus probables.
Construction de la vraisemblance¶
Supposons que nous avons un modèle paramétrique qui, pour chaque entrée et choix de paramètres , définit une distribution sur les sorties possibles . Par exemple, en régression, ce pourrait être une gaussienne centrée sur .
Considérons un seul exemple . Pour des paramètres fixés, nous pouvons évaluer : la probabilité (ou densité) que le modèle assigne à l’observation . Si cette valeur est élevée, les paramètres “expliquent bien” cette observation. Si elle est faible, est une valeur improbable sous ce modèle.
Avec deux exemples indépendants et , la probabilité conjointe est le produit:
Avec exemples indépendants, nous obtenons la vraisemblance:
Cette quantité est une fonction de . Elle répond à la question: pour ce choix de paramètres, quelle est la probabilité d’avoir observé exactement ces données?
Pourquoi maximiser?¶
Si , alors les données observées sont plus probables sous que sous . Les paramètres rendent les observations moins “surprenantes”.
L’estimateur du maximum de vraisemblance (EMV, ou MLE pour maximum likelihood estimator en anglais) choisit les paramètres qui maximisent cette probabilité:
C’est le choix de paramètres sous lequel nos données sont les plus “attendues”.
Du produit à la somme¶
En pratique, multiplier probabilités (souvent petites) pose des problèmes numériques: le résultat devient rapidement trop petit pour être représenté par un ordinateur. Le logarithme résout ce problème: il transforme le produit en somme et, comme c’est une fonction croissante, il ne change pas le maximiseur:
Pour l’optimisation, nous préférons minimiser plutôt que maximiser (par convention). La log-vraisemblance négative (NLV, ou NLL pour negative log-likelihood en anglais) est notre fonction objectif:
Remarquez la structure: c’est une somme sur les exemples d’une quantité qui dépend de chaque observation. Cette quantité joue le rôle d’une fonction de perte. Le maximum de vraisemblance est donc un cas particulier de la minimisation du risque empirique, où la perte est définie par le modèle probabiliste lui-même.
Régression avec bruit gaussien: d’où vient la perte quadratique?¶
Appliquons ce principe à la régression. Le modèle de génération des données suppose que la sortie observée est la prédiction “vraie” du modèle, corrompue par un bruit aléatoire gaussien:
Ce modèle dit que si nous connaissions les vrais paramètres et que nous mesurions pour un donné, nous obtiendrions plus ou moins la plupart du temps.
La distribution conditionnelle qui en découle est:
Calculons la log-vraisemblance négative:
Le second terme ne dépend pas de . Minimiser la NLV revient donc exactement à minimiser la somme des erreurs quadratiques.
C’est un résultat fondamental: la perte quadratique n’est pas un choix arbitraire. Elle découle naturellement de l’hypothèse que les erreurs de mesure suivent une loi gaussienne. Le maximum de vraisemblance sous bruit gaussien coïncide avec les moindres carrés.
Dans ce modèle, nous avons supposé que la variance est constante pour toutes les entrées . C’est ce qu’on appelle la régression homoscédastique (du grec homos, même, et skedasis, dispersion). C’est l’hypothèse standard en régression linéaire.
En pratique, l’incertitude peut varier selon l’entrée. Par exemple, les mesures à haute vitesse peuvent être plus bruitées que celles à basse vitesse. La régression hétéroscédastique modélise cette variation en faisant dépendre la variance de :
où prédit la moyenne et prédit l’écart-type. Ce modèle est plus flexible mais requiert d’apprendre des paramètres supplémentaires.
Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.stats import norm
from IPython.display import HTML
# Générer des données synthétiques
np.random.seed(42)
N = 100
x_data = np.random.uniform(0.5, 9.5, N)
f_mu = lambda x: 0.5 * x + 1
# Homoscédastique: variance constante
sigma_homo = 0.7
y_homo = f_mu(x_data) + np.random.normal(0, sigma_homo, N)
# Hétéroscédastique: variance croissante
f_sigma = lambda x: 0.3 + 0.12 * x
y_hetero = f_mu(x_data) + np.random.normal(0, f_sigma(x_data))
# Configuration de la figure
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
x_line = np.linspace(0, 10, 100)
y_pdf_range = np.linspace(-3, 9, 200)
scale = 2.5 # échelle pour afficher les PDFs
def init():
for ax in axes:
ax.clear()
return []
def animate(frame):
x_current = 0.5 + frame * 9 / 59 # balayer de 0.5 à 9.5
for idx, (ax, y_data, title, color, get_sigma) in enumerate([
(axes[0], y_homo, 'Régression homoscédastique', 'steelblue', lambda x: sigma_homo),
(axes[1], y_hetero, 'Régression hétéroscédastique', 'coral', f_sigma)
]):
ax.clear()
# Données et ligne de régression
ax.scatter(x_data, y_data, alpha=0.4, s=20, c='gray', zorder=1)
ax.plot(x_line, f_mu(x_line), 'k-', linewidth=2, zorder=2)
# Gaussienne à la position actuelle
mu = f_mu(x_current)
sigma = get_sigma(x_current)
pdf = norm.pdf(y_pdf_range, mu, sigma)
# Afficher la gaussienne "horizontalement"
ax.fill_betweenx(y_pdf_range, x_current, x_current + scale * pdf,
alpha=0.5, color=color, zorder=3)
ax.plot(x_current + scale * pdf, y_pdf_range, color=color, linewidth=2, zorder=4)
# Ligne verticale indiquant la position
ax.axvline(x_current, color=color, linestyle='--', alpha=0.5, linewidth=1)
# Point sur la courbe de régression
ax.scatter([x_current], [mu], color='black', s=50, zorder=5)
# Bande ±2σ
ax.fill_between([x_current - 0.1, x_current + 0.1],
[mu - 2*sigma, mu - 2*sigma],
[mu + 2*sigma, mu + 2*sigma],
alpha=0.2, color=color, zorder=0)
ax.set_xlim(-0.5, 12)
ax.set_ylim(-2, 8)
ax.set_xlabel(r'$x$', fontsize=11)
ax.set_ylabel(r'$y$', fontsize=11)
sigma_label = r'$\sigma^2$ constant' if idx == 0 else r'$\sigma^2(x)$ variable'
ax.set_title(f'{title}\n{sigma_label}', fontsize=11)
fig.tight_layout()
return []
anim = FuncAnimation(fig, animate, init_func=init, frames=60, interval=80, blit=True)
anim.save('_static/regression_scedasticity.gif', writer='pillow', fps=12, dpi=100)
plt.close()
# Afficher le GIF
from IPython.display import Image
Image(filename='_static/regression_scedasticity.gif')
L’animation illustre la différence fondamentale entre les deux modèles. À chaque position , la distribution conditionnelle est une gaussienne (la “cloche” colorée) centrée sur la courbe de régression . Dans le cas homoscédastique (gauche), la cloche garde la même largeur partout. Dans le cas hétéroscédastique (droite), la largeur varie avec . Ici, l’incertitude augmente vers la droite, ce qui se traduit par une dispersion plus grande des points.
Régression linéaire homoscédastique¶
Appliquons maintenant le maximum de vraisemblance au cas le plus courant: la régression linéaire homoscédastique. Le modèle probabiliste est:
La fonction de moyenne est linéaire: , où nous avons absorbé le biais en posant . Le vecteur contient donc le biais et les poids . La variance est constante.
Formulation matricielle¶
Avec observations , nous pouvons écrire le modèle sous forme matricielle. Définissons la matrice de conception (design matrix) dont chaque ligne contient une observation augmentée d’un 1 pour le biais:
Le vecteur des prédictions est , et le vecteur des résidus est .
La somme des carrés des résidus (residual sum of squares, RSS) s’écrit:
Comme nous l’avons vu, minimiser la NLV sous bruit gaussien homoscédastique revient à minimiser le RSS (ou de manière équivalente, l’EQM = RSS/).
Solution analytique: les équations normales¶
Pour trouver le minimum, nous calculons le gradient du RSS par rapport à et l’égalons à zéro. En développant:
Le gradient est:
En posant ce gradient égal à zéro, nous obtenons les équations normales:
Si la matrice est inversible, la solution unique est:
Cette solution porte le nom d’estimateur des moindres carrés ordinaires (MCO, ou ordinary least squares, OLS). Elle est exactement équivalente à l’EMV sous l’hypothèse de bruit gaussien homoscédastique.
Conditions d’existence de la solution¶
La matrice est de taille . Elle est inversible si et seulement si est de rang plein colonne, c’est-à-dire si les colonnes de sont linéairement indépendantes. Cela requiert:
Plus d’observations que de paramètres:
Pas de colinéarité parfaite: aucune caractéristique ne doit être une combinaison linéaire exacte des autres
Quand est mal conditionnée (presque singulière), de petites perturbations dans les données peuvent causer de grandes variations dans la solution. C’est l’un des problèmes que la régularisation permet de résoudre.
Source
import numpy as np
import matplotlib.pyplot as plt
# Données synthétiques: relation linéaire avec bruit
np.random.seed(42)
N = 30
x = np.random.uniform(0, 10, N)
y_true = 2.5 * x + 3 # vraie relation: y = 2.5x + 3
y = y_true + np.random.normal(0, 3, N) # ajout de bruit gaussien
# Construction de la matrice de conception (avec colonne de 1 pour le biais)
X = np.column_stack([np.ones(N), x])
# Solution OLS: w = (X^T X)^{-1} X^T y
XtX = X.T @ X
Xty = X.T @ y
w_ols = np.linalg.solve(XtX, Xty)
b_ols, slope_ols = w_ols[0], w_ols[1]
# Prédictions
x_grid = np.linspace(0, 10, 100)
y_pred = b_ols + slope_ols * x_grid
# Visualisation
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
# Gauche: données et droite ajustée
ax = axes[0]
ax.scatter(x, y, alpha=0.7, label='Observations')
ax.plot(x_grid, y_pred, 'k-', linewidth=2, label=f'MCO: $y = {slope_ols:.2f}x + {b_ols:.2f}$')
ax.plot(x_grid, 2.5 * x_grid + 3, 'g--', alpha=0.5, label='Vraie relation')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.legend()
ax.set_title('Régression linéaire par moindres carrés')
# Droite: résidus
ax = axes[1]
y_fitted = b_ols + slope_ols * x
residuals = y - y_fitted
ax.stem(x, residuals, basefmt=' ', linefmt='C0-', markerfmt='C0o')
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel('$x$')
ax.set_ylabel(r'Résidu $r_i = y_i - \hat{y}_i$')
ax.set_title(f'Résidus (RSS = {np.sum(residuals**2):.1f})')
plt.tight_layout()
La figure de gauche montre les données et la droite ajustée par moindres carrés. La figure de droite montre les résidus: les écarts entre les observations et les prédictions. L’EMV minimise la somme des carrés de ces résidus.
Exemple: pharmacocinétique¶
L’EMV s’applique à des modèles non linéaires. Considérons la concentration d’un médicament dans le sang après administration orale. Les données suivantes proviennent d’une étude sur la théophylline, un bronchodilatateur:
Source
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# Données pharmacocinétiques: théophylline, sujet 1 (Boeckmann et al., 1994)
time = np.array([0, 0.25, 0.57, 1.12, 2.02, 3.82, 5.10, 7.03, 9.05, 12.12, 24.37])
conc = np.array([0.74, 2.84, 6.57, 10.50, 9.66, 8.58, 8.36, 7.47, 6.89, 5.94, 3.28])
# Model: C(t) = C0 * exp(-k * t) for t > t_peak
# We'll fit on the decay phase (after peak)
peak_idx = np.argmax(conc)
t_decay = time[peak_idx:]
c_decay = conc[peak_idx:]
# MLE: minimize NLL under Gaussian noise
def neg_log_likelihood(params, t, c):
C0, k, sigma = params
if sigma <= 0 or k <= 0 or C0 <= 0:
return np.inf
pred = C0 * np.exp(-k * (t - t[0]))
nll = 0.5 * len(t) * np.log(2 * np.pi * sigma**2)
nll += 0.5 * np.sum((c - pred)**2) / sigma**2
return nll
# Initial guess and optimization
x0 = [c_decay[0], 0.1, 1.0]
result = minimize(neg_log_likelihood, x0, args=(t_decay, c_decay), method='Nelder-Mead')
C0_mle, k_mle, sigma_mle = result.x
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
# Left: data and fit
ax = axes[0]
ax.scatter(time, conc, s=50, zorder=5, label='Observations')
t_grid = np.linspace(time[peak_idx], 25, 100)
ax.plot(t_grid, C0_mle * np.exp(-k_mle * (t_grid - time[peak_idx])), 'k--',
label=f'EMV: $C_0$={C0_mle:.1f}, $k$={k_mle:.2f}')
ax.axvline(time[peak_idx], color='gray', linestyle=':', alpha=0.5)
ax.set_xlabel('Temps (h)')
ax.set_ylabel('Concentration (mg/L)')
ax.legend()
ax.set_title('Concentration plasmatique de théophylline')
# Right: residuals
ax = axes[1]
pred_decay = C0_mle * np.exp(-k_mle * (t_decay - t_decay[0]))
residuals = c_decay - pred_decay
ax.stem(t_decay, residuals, basefmt=' ')
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel('Temps (h)')
ax.set_ylabel('Résidu (mg/L)')
ax.set_title(rf'$\sigma$ estimé: {sigma_mle:.2f} mg/L')
plt.tight_layout()
Le modèle décrit la décroissance exponentielle après le pic de concentration. Les paramètres (concentration initiale) et (constante d’élimination) sont estimés par maximum de vraisemblance sous l’hypothèse d’un bruit gaussien. Cette approche est identique à celle des moindres carrés, mais elle fournit également une estimation de l’écart-type du bruit .
Classification binaire¶
La perte 0-1 pour la classification est discontinue, ce qui empêche l’utilisation de méthodes de gradient. La fonction sigmoïde contourne ce problème: c’est une approximation lisse de la fonction échelon (step function). Elle transforme n’importe quel score réel en une valeur dans l’intervalle , que nous pouvons interpréter comme une probabilité.
Source
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# Create figure
fig, ax = plt.subplots(figsize=(8, 5))
# Define x range
z = np.linspace(-4, 4, 200)
# Step function (Heaviside)
step = (z >= 0).astype(float)
# Sigmoid function with temperature parameter
def sigmoid(z, alpha=1):
return 1 / (1 + np.exp(-alpha * z))
# Initialize plot
line_step, = ax.plot(z, step, 'k--', linewidth=2, label='Fonction échelon', alpha=0.7)
line_sigmoid, = ax.plot([], [], 'b-', linewidth=2, label='Sigmoïde $\\sigma(\\alpha z)$')
ax.axhline(0.5, color='gray', linestyle=':', alpha=0.5, linewidth=1)
ax.axvline(0, color='gray', linestyle=':', alpha=0.5, linewidth=1)
ax.set_xlim(-4, 4)
ax.set_ylim(-0.1, 1.1)
ax.set_xlabel('$z$')
ax.set_ylabel('$\\sigma(\\alpha z)$')
ax.set_title('Approximation de la fonction échelon par la sigmoïde')
ax.legend(loc='best')
ax.grid(True, alpha=0.3)
# Animation function
def animate(frame):
# Alpha increases from 0.5 to 10
alpha = 0.5 + (frame / 100) * 9.5
y = sigmoid(z, alpha)
line_sigmoid.set_data(z, y)
ax.set_title(f'Approximation de la fonction échelon par la sigmoïde ($\\alpha = {alpha:.2f}$)')
return line_sigmoid,
# Create animation
anim = FuncAnimation(fig, animate, frames=100, interval=50, blit=True, repeat=True)
anim.save('_static/sigmoid_approximation.gif', writer='pillow', fps=20, dpi=100)
plt.close()
# Afficher le GIF
from IPython.display import Image
Image(filename='_static/sigmoid_approximation.gif')
L’animation montre comment la sigmoïde se rapproche de la fonction échelon lorsque le paramètre augmente. Pour , la sigmoïde est douce; pour grand, elle devient presque aussi abrupte que la fonction échelon, tout en restant différentiable.
Cette interprétation probabiliste n’est pas qu’une astuce numérique. Elle correspond exactement à modéliser par une distribution de Bernoulli dont le paramètre dépend de l’entrée.
Pour la classification binaire avec , nous modélisons la probabilité de la classe positive par:
où est la fonction sigmoïde et est le logit (ou log-odds), le score brut du modèle avant transformation. Le logit est le logarithme du rapport des probabilités: . La distribution conditionnelle suit une loi de Bernoulli:
La log-vraisemblance négative est:
Cette quantité est l’entropie croisée binaire. Elle correspond à la perte logistique, à une reparamétrisation près.
Classification multiclasse¶
Pour la classification avec classes (), nous généralisons le modèle binaire en utilisant la distribution catégorielle (ou multinomiale). Au lieu de modéliser une seule probabilité , nous modélisons un vecteur de probabilités où et .
Pour transformer les scores bruts du modèle en probabilités, nous utilisons la fonction softmax:
où est le score pour la classe . La fonction softmax généralise la sigmoïde au cas multiclasse: elle transforme scores réels en un vecteur de probabilités qui somme à 1.
La distribution conditionnelle suit une loi catégorielle:
où vaut 1 si et 0 sinon. En utilisant l’encodage one-hot , cette expression devient:
La log-vraisemblance négative est:
où . Cette quantité est l’entropie croisée multiclasse. Elle généralise l’entropie croisée binaire au cas où il y a plus de deux classes.
Pour la classification binaire avec , le softmax se réduit à la sigmoïde. En effet, si nous définissons , alors:
Le modèle binaire et le modèle multiclasse partagent donc la même structure probabiliste, avec la distribution catégorielle comme généralisation naturelle de la distribution de Bernoulli.
Maximum a posteriori¶
Plutôt que de travailler avec la distribution a posteriori complète (ce qui peut être coûteux), nous pouvons chercher son mode: la valeur des paramètres la plus probable a posteriori. C’est l’estimateur du maximum a posteriori (MAP):
Le dénominateur ne dépend pas de et peut être ignoré pour l’optimisation. En passant au logarithme:
Cette expression révèle une structure familière. Si nous posons , nous obtenons:
C’est exactement la forme du risque empirique régularisé. La régularisation correspond à l’ajout d’un a priori sur les paramètres. Le terme de régularisation est le logarithme négatif de la distribution a priori.
Le maximum de vraisemblance comme cas particulier¶
Que se passe-t-il si nous n’avons aucune préférence a priori sur les paramètres? Cela correspond à un a priori uniforme (ou constant): .
Dans ce cas, est une constante qui n’affecte pas l’optimisation, et le MAP se réduit à l’EMV:
L’EMV est donc un cas particulier du MAP: celui où nous supposons implicitement que toutes les valeurs de paramètres sont également plausibles avant d’observer les données. Cette perspective unifie les deux approches dans un même cadre.
Limites de l’a priori uniforme¶
L’a priori uniforme (et donc l’EMV) peut être problématique quand les données sont peu nombreuses. Considérons l’estimation de la probabilité qu’une pièce tombe sur face.
Supposons que nous lancions la pièce 3 fois et obtenions 3 faces. L’estimateur du maximum de vraisemblance pour une distribution de Bernoulli est:
où est le nombre de faces et le nombre de piles. Cette estimation dit que la probabilité d’obtenir face est de 100%. Si nous utilisions ce modèle pour prédire de futurs lancers, nous prédirons toujours face, ce qui est peu plausible pour une vraie pièce.
Le problème est que l’EMV (avec son a priori uniforme implicite) dispose de suffisamment de flexibilité pour reproduire parfaitement les données d’entraînement, même quand celles-ci sont peu nombreuses ou non représentatives. Un a priori informatif peut atténuer ce problème.
Source
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# MLE vs MAP for Bernoulli with few observations
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
# Different sample sizes
samples_list = [
[1, 1, 1], # 3 heads
[1, 1, 1, 0], # 3 heads, 1 tail
[1, 1, 1, 0, 0, 1, 1, 0, 1, 1] # 7 heads, 3 tails
]
theta_grid = np.linspace(0.001, 0.999, 200)
for ax, samples in zip(axes, samples_list):
n1 = sum(samples) # heads
n0 = len(samples) - n1 # tails
# MLE
theta_mle = n1 / (n0 + n1)
# Likelihood (unnormalized)
likelihood = theta_grid**n1 * (1 - theta_grid)**n0
likelihood = likelihood / likelihood.max()
ax.plot(theta_grid, likelihood, 'b-', linewidth=2, label='Vraisemblance')
ax.axvline(theta_mle, color='b', linestyle='--', alpha=0.7,
label=f'EMV: {theta_mle:.2f}')
ax.axvline(0.5, color='gray', linestyle=':', alpha=0.5, label=r'$\theta = 0.5$')
ax.set_xlabel(r'$\theta$')
ax.set_ylabel('Vraisemblance (normalisée)')
ax.set_title(f'{n1} faces, {n0} piles (N={len(samples)})')
ax.legend(fontsize=8)
ax.set_xlim(0, 1)
plt.tight_layout()
La figure montre la vraisemblance pour différents échantillons. Avec seulement 3 observations (toutes faces), la vraisemblance est maximale à . En augmentant la taille de l’échantillon, l’estimation devient plus raisonnable. Voyons comment un a priori non uniforme peut aider.
Exemple: lissage de Laplace¶
Revenons à notre exemple de la pièce de monnaie. Utilisons un a priori Beta sur :
Les paramètres et contrôlent la forme de l’a priori. Pour , l’a priori favorise des valeurs de proches de 0.5.
Le logarithme de l’a posteriori (vraisemblance plus a priori) est:
En dérivant et en résolvant, l’estimateur MAP est:
Avec et nos 3 observations de faces:
Cette estimation est plus raisonnable que l’EMV . L’a priori “tire” l’estimation vers des valeurs moins extrêmes.
Le choix correspond au lissage de Laplace (ou add-one smoothing): c’est comme si nous avions observé une face et une pile supplémentaires avant de commencer. Cette technique est particulièrement utile quand certains événements n’ont jamais été observés dans les données.
Source
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
theta_grid = np.linspace(0.001, 0.999, 200)
# Left: Prior, likelihood, posterior
ax = axes[0]
n1, n0 = 3, 0 # 3 heads, 0 tails
a, b = 2, 2 # Beta prior parameters
# Prior
prior = stats.beta.pdf(theta_grid, a, b)
prior = prior / prior.max()
# Likelihood
likelihood = theta_grid**n1 * (1 - theta_grid)**n0
likelihood = likelihood / likelihood.max()
# Posterior (Beta(a + n1, b + n0))
posterior = stats.beta.pdf(theta_grid, a + n1, b + n0)
posterior = posterior / posterior.max()
ax.plot(theta_grid, prior, 'g-', linewidth=2, label='A priori Beta(2,2)')
ax.plot(theta_grid, likelihood, 'b--', linewidth=2, label='Vraisemblance')
ax.plot(theta_grid, posterior, 'r-', linewidth=2, label='A posteriori')
theta_mle = n1 / (n1 + n0)
theta_map = (n1 + a - 1) / (n1 + n0 + a + b - 2)
ax.axvline(theta_mle, color='b', linestyle=':', alpha=0.7, label=f'EMV: {theta_mle:.2f}')
ax.axvline(theta_map, color='r', linestyle=':', alpha=0.7, label=f'MAP: {theta_map:.2f}')
ax.set_xlabel(r'$\theta$')
ax.set_ylabel('Densité (normalisée)')
ax.set_title('3 faces, 0 pile')
ax.legend(fontsize=8)
ax.set_xlim(0, 1)
# Right: Effect of different priors
ax = axes[1]
priors = [(1, 1, 'Uniforme'), (2, 2, 'Beta(2,2)'), (5, 5, 'Beta(5,5)')]
for a, b, label in priors:
theta_map = (n1 + a - 1) / (n1 + n0 + a + b - 2)
posterior = stats.beta.pdf(theta_grid, a + n1, b + n0)
posterior = posterior / posterior.max()
ax.plot(theta_grid, posterior, linewidth=2, label=f'{label}: MAP={theta_map:.2f}')
ax.axvline(1.0, color='gray', linestyle='--', alpha=0.5, label='EMV: 1.00')
ax.set_xlabel(r'$\theta$')
ax.set_ylabel('A posteriori (normalisé)')
ax.set_title('Effet de différents a priori')
ax.legend(fontsize=8)
ax.set_xlim(0, 1)
plt.tight_layout()
La figure de gauche montre comment l’a posteriori combine l’a priori et la vraisemblance. L’a priori Beta(2,2) “tire” l’estimation vers 0.5, résultant en un MAP de 0.8 au lieu de l’EMV de 1.0. La figure de droite montre l’effet de différents a priori: plus l’a priori est fort (variance faible), plus l’estimation est proche de 0.5.
Régression ridge = MAP avec a priori gaussien¶
Appliquons maintenant ce cadre bayésien à la régression linéaire. Si nous plaçons un a priori gaussien isotrope sur les paramètres:
cet a priori exprime la croyance que les paramètres sont probablement proches de zéro, avec une incertitude contrôlée par .
Le logarithme négatif de cet a priori est:
L’estimateur MAP devient:
C’est exactement la régression ridge, avec . Cette correspondance nous donne une interprétation de l’hyperparamètre:
Grande valeur de (petite variance ): forte croyance que les paramètres sont proches de zéro
Petite valeur de (grande variance ): a priori peu informatif, on fait confiance aux données
L’a priori gaussien sur les paramètres est parfois appelé dégradation des poids (weight decay) dans le contexte des réseaux de neurones, car il “tire” les paramètres vers zéro pendant l’entraînement.
Unification: deux langages pour un même problème¶
Les sections précédentes ont présenté deux approches pour l’apprentissage supervisé. La première, fondée sur la théorie de la décision, définit une fonction de perte et minimise le risque empirique. La seconde, probabiliste, modélise la distribution des données et estime les paramètres par maximum de vraisemblance ou maximum a posteriori.
Ces deux approches semblent différentes, mais elles aboutissent aux mêmes algorithmes. En choisissant la perte logarithmique , le risque empirique devient exactement la log-vraisemblance négative (à un facteur près). Minimiser l’un revient à minimiser l’autre. Sous bruit gaussien, cette perte se réduit à la perte quadratique; sous modèle de Bernoulli, à l’entropie croisée.
De même, ajouter une régularisation au risque empirique revient à supposer un a priori gaussien sur les paramètres. La régression ridge n’est rien d’autre que l’estimation MAP avec cet a priori. Le coefficient encode la force de notre croyance a priori: plus est grand, plus nous “tirons” les paramètres vers zéro.
Pourquoi alors utiliser deux langages? Parce qu’ils éclairent des aspects différents du problème. Le langage décisionnel (risque, perte, minimisation) est opérationnel: il dit comment construire un algorithme. Le langage probabiliste (vraisemblance, a priori, a posteriori) est interprétatif: il dit ce que nous supposons sur les données et pourquoi nos choix sont raisonnables. Ensemble, ils permettent de concevoir des algorithmes et de comprendre leur comportement.
Interprétation informationnelle¶
La théorie de l’information offre une troisième perspective. L’EMV peut se comprendre comme la recherche du modèle paramétrique le plus proche de la distribution empirique des données.
La distribution empirique place une masse sur chaque observation:
La divergence de Kullback-Leibler mesure la dissimilarité entre deux distributions:
Cette quantité est toujours positive ou nulle, et vaut zéro si et seulement si les deux distributions sont identiques. En posant (ce que nous avons observé) et (notre modèle), on peut montrer que:
L’EMV trouve les paramètres qui rendent notre modèle aussi proche que possible de ce que nous avons observé, au sens de la divergence KL. Cette interprétation géométrique complète les perspectives décisionnelle et probabiliste.
Résumé¶
Ce chapitre a développé le cadre formel de l’apprentissage supervisé à travers deux perspectives complémentaires.
La perspective décisionnelle définit le risque comme l’erreur moyenne attendue sur de nouvelles données, puis le remplace par le risque empirique (calculable sur les données d’entraînement). La minimisation du risque empirique cherche le modèle qui minimise cette erreur d’entraînement, en espérant qu’il généralisera. La régularisation pénalise la complexité pour éviter le surapprentissage.
La perspective probabiliste modélise explicitement comment les données ont été générées. Le maximum de vraisemblance (EMV) trouve les paramètres qui rendent les observations les plus probables; le maximum a posteriori (MAP) incorpore des croyances a priori sur les paramètres.
Ces deux perspectives aboutissent aux mêmes algorithmes. Sous bruit gaussien, l’EMV coïncide avec les moindres carrés. La régularisation correspond au MAP avec a priori gaussien sur les paramètres. Le langage décisionnel est opérationnel (comment construire l’algorithme); le langage probabiliste est interprétatif (pourquoi ces choix sont raisonnables).
Le chapitre suivant développe les outils théoriques pour quantifier quand et comment le risque empirique prédit le vrai risque.
Exercices¶
Exercice 1: Usure d’outil
Un machiniste mesure l’usure d’un outil de coupe (en mm) à différents temps de coupe (en minutes):
import numpy as np
# Données d'usure d'outil (simulées selon une loi de puissance avec bruit)
time = np.array([2.0, 5.1, 8.2, 11.3, 14.4, 17.6, 20.7, 23.8, 26.9, 30.0])
wear = np.array([0.08, 0.14, 0.17, 0.21, 0.22, 0.25, 0.27, 0.28, 0.31, 0.32])L’outil doit être remplacé lorsque l’usure atteint 0.4 mm.
Visualisation. Tracez les données. Quelle forme de relation observez-vous?
Ajustement. Ajustez un modèle linéaire et un modèle en loi de puissance aux données. Pour le second modèle, utilisez une transformation logarithmique: .
Comparaison. Calculez l’EQM de chaque modèle sur les données. Lequel ajuste mieux?
Prédiction. Selon chaque modèle, à quel moment l’usure atteindra-t-elle 0.4 mm? Les deux modèles donnent-ils la même réponse?
Extrapolation. Si vous n’aviez mesuré que jusqu’à min, vos prédictions changeraient-elles? Discutez du risque d’extrapolation.
Solution Exercice 1
Visualisation. Les données montrent une relation non linéaire, concave: l’usure augmente rapidement au début puis ralentit. Cela suggère une loi de puissance avec exposant .
Ajustement.
Modèle linéaire:
coeffs = np.polyfit(time, wear, 1)donne , .Loi de puissance: en posant , on ajuste une droite dans l’espace log-log:
coeffs = np.polyfit(np.log(time), np.log(wear), 1). On obtient et .
Comparaison. L’EQM du modèle linéaire est typiquement plus élevé car il ne capture pas la courbure. Le modèle en loi de puissance ajuste mieux les données.
Prédiction. Pour trouver tel que :
Linéaire:
Puissance:
Les réponses diffèrent significativement car les modèles extrapolent différemment.
Extrapolation. Avec moins de données, les estimations des paramètres changent, et les prédictions au-delà des données observées deviennent plus incertaines. L’extrapolation est risquée car le comportement futur peut ne pas suivre le modèle ajusté sur les données passées.
Exercice 2: Risque et risque empirique
Soit un problème de classification binaire avec la perte 0-1. Un classificateur fait 3 erreurs sur 20 exemples d’entraînement.
Quel est le risque empirique de sur l’ensemble d’entraînement?
Peut-on en déduire le vrai risque ? Pourquoi ou pourquoi pas?
Si nous avions 1000 exemples de test et que fait 45 erreurs, quelle serait notre meilleure estimation du vrai risque?
Solution Exercice 2
Risque empirique: (soit 15% d’erreur).
Non, on ne peut pas en déduire le vrai risque. Le risque empirique sur l’entraînement est une estimation biaisée du vrai risque car:
Le modèle a été choisi/optimisé pour bien performer sur ces mêmes données
Il y a surapprentissage potentiel: peut avoir mémorisé des particularités de l’entraînement qui ne généralisent pas
Le risque empirique sur l’entraînement sous-estime généralement le vrai risque
Estimation sur le test: (soit 4.5% d’erreur). Cette estimation est plus fiable car:
Les données de test n’ont pas été utilisées pour construire
Avec 1000 exemples, l’estimation est plus précise (écart-type )
Exercice 3: Maximum de vraisemblance
Soit un échantillon i.i.d. d’une distribution exponentielle de paramètre :
Écrivez la vraisemblance et la log-vraisemblance .
Dérivez l’estimateur du maximum de vraisemblance (EMV) .
Si les observations sont , calculez .
Solution Exercice 3
Exercice 4: Fonctions de perte
Soit (classe positive) et un score .
Calculez la perte 0-1, la perte logistique, et la perte à charnière.
Répétez pour (prédiction incorrecte).
Tracez les trois fonctions de perte en fonction de pour . Vérifiez que les pertes de substitution majorent la perte 0-1.
Solution Exercice 4
Pour et (prédiction correcte, marge ):
Perte 0-1:
Perte logistique:
Perte à charnière:
Pour et (prédiction incorrecte, marge ):
Perte 0-1:
Perte logistique:
Perte à charnière:
Vérification graphique:
import numpy as np import matplotlib.pyplot as plt margin = np.linspace(-3, 3, 100) loss_01 = (margin < 0).astype(float) loss_log = np.log(1 + np.exp(-margin)) loss_hinge = np.maximum(0, 1 - margin) plt.plot(margin, loss_01, label='0-1') plt.plot(margin, loss_log, label='Logistique') plt.plot(margin, loss_hinge, label='Charnière') plt.legend()On vérifie que pour tout : et .
Exercice 5: Dérivation des moindres carrés ordinaires ★
Soit un problème de régression linéaire simple avec observations:
Écrivez la somme des carrés des résidus .
Calculez les dérivées partielles et .
En posant ces dérivées égales à zéro, résolvez le système d’équations pour obtenir les estimateurs et .
Application numérique: Pour les données , calculez les coefficients MCO et , puis la prédiction pour .
Solution Exercice 5
Exercice 6: Maximum de vraisemblance pour classification binaire ★
Soit un problème de classification binaire avec observations où . On modélise la probabilité de la classe positive par:
Écrivez la vraisemblance pour observations i.i.d. suivant une distribution de Bernoulli.
Écrivez la log-vraisemblance .
Montrez que maximiser la log-vraisemblance revient à minimiser l’entropie croisée binaire.
Application: Pour les observations avec les scores égaux à , calculez la log-vraisemblance.
Solution Exercice 6
Exercice 7: Expansion de caractéristiques ★
Considérez un problème de régression où la relation entre et est quadratique:
Générez points selon ce modèle pour .
Ajustez un modèle linéaire aux données. Visualisez le résultat et calculez l’EQM.
Définissez la transformation . Montrez que le problème devient linéaire en dans cet espace étendu.
Ajustez un modèle linéaire dans l’espace étendu: . Comparez l’EQM avec le modèle précédent.
Que se passe-t-il si vous utilisez ? Discutez du risque de surapprentissage.
Solution Exercice 7
Génération des données:
import numpy as np np.random.seed(42) x = np.random.uniform(-2, 4, 50) y = 3 - 2*x + 0.5*x**2 + np.random.normal(0, 0.5, 50)Modèle linéaire simple:
coeffs_lin = np.polyfit(x, y, 1) y_pred_lin = np.polyval(coeffs_lin, x) mse_lin = np.mean((y - y_pred_lin)**2)L’EQM est élevé car le modèle linéaire ne peut pas capturer la courbure. Visuellement, la droite ne suit pas la tendance parabolique des données.
Transformation en problème linéaire:
En définissant , le modèle s’écrit:
C’est linéaire en , même si c’est non-linéaire en . On peut appliquer les équations normales dans cet espace.
Modèle quadratique:
coeffs_quad = np.polyfit(x, y, 2) y_pred_quad = np.polyval(coeffs_quad, x) mse_quad = np.mean((y - y_pred_quad)**2)L’EQM est beaucoup plus faible (proche de ). Les coefficients récupérés sont proches des vrais: , , .
Expansion de degré 5:
Avec plus de termes, le modèle peut s’ajuster encore mieux aux données d’entraînement (EQM très faible), mais il risque de capturer le bruit plutôt que le signal. Sur de nouvelles données, les performances se dégradent. C’est le surapprentissage: le modèle a trop de capacité par rapport à la complexité réelle de la relation.
Exercice 8: EMV comme cas particulier de MAP ★
Le théorème de Bayes nous donne la distribution a posteriori des paramètres:
L’estimateur MAP maximise cette distribution a posteriori.
Écrivez la log-posterior en fonction de la log-vraisemblance et du log a priori.
Supposons un a priori uniforme (constant): pour tout . Montrez que l’estimateur MAP se réduit à l’estimateur du maximum de vraisemblance (EMV).
Pour quels autres types d’a priori l’EMV et le MAP coïncident-ils?
Expliquez pourquoi utiliser un a priori uniforme peut être problématique dans certains cas.
Solution Exercice 8
Log-posterior:
Le terme ne dépend pas de , donc pour l’optimisation:
A priori uniforme:
Si (constante), alors est aussi une constante.
L’EMV est donc un cas particulier du MAP avec a priori uniforme.
Autres a priori:
L’EMV et le MAP coïncident pour tout a priori qui est constant sur le domaine des paramètres, ou plus généralement, pour tout a priori dont le log est constant (à une constante additive près). Cela inclut les a priori impropres (non normalisables) qui sont uniformes sur .
Problèmes de l’a priori uniforme:
Avec peu de données: l’EMV peut être extrême. Exemple: 3 lancers de pièce donnant 3 faces → EMV = 100% de probabilité de face.
Paramètres non bornés: un a priori uniforme sur n’est pas une vraie distribution de probabilité (a priori impropre).
Invariance: un a priori uniforme sur n’est pas uniforme sur pour une transformation non-linéaire .
Pas d’information: on ignore toute connaissance préalable sur les valeurs plausibles des paramètres.
Exercice 9: Régression ridge et colinéarité ★★
La colinéarité entre les caractéristiques rend la matrice mal conditionnée, ce qui peut déstabiliser la solution MCO.
Générez des données avec deux caractéristiques presque colinéaires:
np.random.seed(42) n = 30 x1 = np.random.randn(n) x2 = x1 + 0.01 * np.random.randn(n) # x2 ≈ x1 y = 2*x1 + 3*x2 + 0.5*np.random.randn(n)Calculez le nombre de conditionnement de (avec
np.linalg.cond). Que signifie un grand nombre de conditionnement?Ajustez un modèle MCO. Les coefficients et sont-ils proches des vraies valeurs (2 et 3)?
Ajustez des modèles Ridge pour . Comment les coefficients évoluent-ils?
Tracez le “chemin de régularisation”: les coefficients en fonction de .
Solution Exercice 9
Génération des données: (code fourni dans l’énoncé)
Nombre de conditionnement:
X = np.column_stack([np.ones(n), x1, x2]) cond = np.linalg.cond(X.T @ X) print(f"Conditionnement: {cond:.0f}")Le nombre de conditionnement est très élevé (de l’ordre de 106 ou plus). Cela signifie que de petites perturbations dans les données peuvent causer de grandes variations dans la solution. La matrice est proche d’être singulière.
Modèle MCO:
theta_ols = np.linalg.solve(X.T @ X, X.T @ y)Les coefficients MCO sont très instables: et peuvent être très différents de 2 et 3, et parfois de signes opposés avec de grandes magnitudes. Le modèle “distribue” l’effet entre les deux variables de manière arbitraire.
Modèles Ridge:
from sklearn.linear_model import Ridge for lam in [0.01, 0.1, 1, 10]: model = Ridge(alpha=lam, fit_intercept=True) model.fit(np.column_stack([x1, x2]), y) print(f"λ={lam}: θ1={model.coef_[0]:.2f}, θ2={model.coef_[1]:.2f}")Avec croissant, les coefficients se rapprochent de zéro et deviennent plus stables. Les deux coefficients convergent vers des valeurs similaires (autour de 2.5 chacun), ce qui reflète mieux la symétrie du problème.
Chemin de régularisation:
lambdas = np.logspace(-3, 2, 50) coefs = [] for lam in lambdas: model = Ridge(alpha=lam, fit_intercept=True) model.fit(np.column_stack([x1, x2]), y) coefs.append(model.coef_) coefs = np.array(coefs) plt.plot(np.log10(lambdas), coefs[:, 0], label='θ1') plt.plot(np.log10(lambdas), coefs[:, 1], label='θ2') plt.xlabel('log10(λ)') plt.ylabel('Coefficients') plt.legend()On observe que pour petit, les coefficients sont instables et peuvent être extrêmes. Pour grand, ils convergent vers zéro. Il existe une zone intermédiaire où les coefficients sont raisonnables.
Exercice 10: Softmax et classification multiclasse ★★
La fonction softmax transforme un vecteur de scores en un vecteur de probabilités:
Pour classes et les scores , calculez manuellement . Vérifiez que les probabilités somment à 1.
Montrez que le softmax est invariant par translation: pour tout scalaire .
Pour classes, montrez que le softmax se réduit à la sigmoïde. Posez et montrez que .
Pour un problème à 3 classes avec les vraies étiquettes one-hot (classe 2) et les probabilités prédites , calculez l’entropie croisée.
Solution Exercice 10
Calcul du softmax:
Vérification: ✓
Invariance par translation:
Cette propriété est utile numériquement: on peut soustraire pour éviter les débordements.
Cas :
En posant , on retrouve bien .
Entropie croisée:
Seule la composante correspondant à la vraie classe contribue à la perte.
Exercice 11: MAP avec a priori gaussien et régression ridge ★★
Considérons un modèle de régression linéaire gaussien:
avec un a priori gaussien isotrope sur les paramètres:
Écrivez la log-vraisemblance pour observations i.i.d.
Écrivez le log a priori .
Montrez que l’estimateur MAP s’écrit:
et identifiez en fonction de et .
Interprétez: que signifie un grand ? Un petit ?
Solution Exercice 11
Log-vraisemblance:
Log a priori:
Estimateur MAP:
En ignorant les constantes:
En multipliant par (qui ne change pas l’argmax):
Donc .
Interprétation:
Grand (a priori large): on est peu sûr que les paramètres sont proches de zéro, donc petit, peu de régularisation, MAP proche de EMV.
Petit (a priori concentré): on croit fortement que les paramètres sont proches de zéro, donc grand, forte régularisation, coefficients tirés vers zéro.
Le rapport compare l’incertitude dans les données () à l’incertitude dans l’a priori (). Plus les données sont bruitées, plus on fait confiance à l’a priori.
Exercice 12: Homoscédasticité et hétéroscédasticité ★★
En régression, l’homoscédasticité suppose que la variance du bruit est constante: . L’hétéroscédasticité suppose que la variance dépend de .
Générez deux jeux de données (, ):
Homoscédastique: avec
Hétéroscédastique: avec
Visualisez les deux jeux de données. Quelle différence observez-vous?
Ajustez un modèle linéaire sur chaque jeu. Les coefficients sont-ils similaires?
Tracez les résidus en fonction de pour les deux cas. Que remarquez-vous?
Pourquoi l’hétéroscédasticité peut-elle être problématique pour l’inférence statistique (intervalles de confiance, tests)?
Solution Exercice 12
Génération des données:
np.random.seed(42) x = np.random.uniform(0, 10, 100) # Homoscédastique y_homo = 2*x + np.random.normal(0, 1, 100) # Hétéroscédastique y_hetero = 2*x + np.random.normal(0, 0.3*x, 100)Visualisation:
Dans le cas homoscédastique, les points sont dispersés uniformément autour de la droite sur toute la plage de . Dans le cas hétéroscédastique, la dispersion augmente avec : les points sont serrés près de et très dispersés pour les grandes valeurs de .
Coefficients:
Les coefficients MCO sont similaires dans les deux cas (proches de , ). MCO reste non biaisé sous hétéroscédasticité, mais n’est plus optimal (pas de variance minimale).
Résidus:
Homoscédastique: les résidus sont répartis uniformément autour de zéro, avec une dispersion constante.
Hétéroscédastique: les résidus montrent un “cône” ou “éventail” (fan shape): la dispersion augmente avec . C’est le signe classique d’hétéroscédasticité.
Problèmes d’inférence:
Les erreurs standard des coefficients sont incorrectes: elles supposent une variance constante.
Les intervalles de confiance et tests t ne sont pas valides.
Les tests de significativité peuvent être trop optimistes ou trop pessimistes.
Solution: utiliser des erreurs standard robustes (Huber-White) ou des moindres carrés pondérés.
Exercice 13: Validation croisée pour le choix de λ ★★
La validation croisée permet de choisir l’hyperparamètre de la régression ridge sans utiliser de données de test.
Générez un jeu de données de régression polynomiale ():
np.random.seed(42) x = np.random.uniform(-3, 3, 50) y = 0.5*x**3 - x**2 + 2*x + np.random.normal(0, 2, 50)Créez une matrice de caractéristiques polynomiales de degré 10: .
Implémentez la validation croisée à 5 plis (5-fold CV):
Divisez les données en 5 groupes
Pour chaque :
Entraînez sur 4 plis, évaluez sur le 5ème
Calculez l’EQM moyen sur les 5 plis
Tracez l’EQM de validation en fonction de . Quel choisiriez-vous?
Comparez les performances (sur un ensemble de test séparé) de MCO () et de Ridge avec le optimal.
Solution Exercice 13
1-2. Génération et matrice de caractéristiques:
from sklearn.preprocessing import PolynomialFeatures
np.random.seed(42)
x = np.random.uniform(-3, 3, 50)
y = 0.5*x**3 - x**2 + 2*x + np.random.normal(0, 2, 50)
poly = PolynomialFeatures(degree=10, include_bias=True)
X = poly.fit_transform(x.reshape(-1, 1))Validation croisée:
from sklearn.model_selection import KFold from sklearn.linear_model import Ridge lambdas = np.logspace(-4, 2, 20) kf = KFold(n_splits=5, shuffle=True, random_state=42) cv_scores = [] for lam in lambdas: fold_scores = [] for train_idx, val_idx in kf.split(X): model = Ridge(alpha=lam, fit_intercept=False) model.fit(X[train_idx], y[train_idx]) y_pred = model.predict(X[val_idx]) mse = np.mean((y[val_idx] - y_pred)**2) fold_scores.append(mse) cv_scores.append(np.mean(fold_scores))Visualisation et choix de λ:
plt.plot(np.log10(lambdas), cv_scores) plt.xlabel('log10(λ)') plt.ylabel('EQM de validation') best_idx = np.argmin(cv_scores) best_lambda = lambdas[best_idx]La courbe montre typiquement:
EQM élevé pour très petit (surapprentissage)
EQM minimal pour intermédiaire
EQM qui remonte pour grand (sous-apprentissage)
Le optimal se situe au minimum de la courbe (souvent autour de 10-1 à 100).
Comparaison finale:
MCO avec degré 10 surapprend fortement et a un EQM de test élevé. Ridge avec optimal a un EQM de test beaucoup plus faible car la régularisation empêche les coefficients d’exploser.
````{admonition} Exercice 14: Prédicteur de Bayes optimal (perte quadratique) ★★
:class: hint dropdown
Le prédicteur de Bayes optimal minimise le risque pour une perte donnée, en supposant que la vraie distribution $p(y|x)$ est connue.
Considérons la distribution conditionnelle suivante pour un $x$ donné:
$$
p(y|x) = 0.3 \cdot \mathcal{N}(y | 1, 0.5^2) + 0.7 \cdot \mathcal{N}(y | 4, 1^2)
$$
C'est un mélange de deux gaussiennes.
1. Tracez cette distribution $p(y|x)$.
2. Calculez l'espérance $\mathbb{E}[y|x]$ (utiliser la linéarité de l'espérance).
3. Pour la perte quadratique, le prédicteur optimal est la moyenne conditionnelle. Quelle est donc la prédiction optimale $\hat{y}^*$?
4. Calculez le risque de Bayes (l'erreur minimale atteignable): $\mathcal{R}^* = \mathbb{E}[(y - \hat{y}^*)^2 | x]$.
5. Si vous prédisiez le mode (la valeur la plus probable) au lieu de la moyenne, quelle serait votre prédiction? Quel serait le risque?Solution Exercice 14
Distribution:
y = np.linspace(-2, 8, 200) p = 0.3 * scipy.stats.norm.pdf(y, 1, 0.5) + 0.7 * scipy.stats.norm.pdf(y, 4, 1) plt.plot(y, p) plt.xlabel('y') plt.ylabel('p(y|x)')La distribution est bimodale avec un petit pic à et un grand pic à .
Espérance:
Par linéarité:
Prédiction optimale:
Pour la perte quadratique,
Risque de Bayes:
Pour un mélange:
où est la composante.
Prédiction par le mode:
Le mode est le pic le plus haut, soit (puisque le poids 0.7 > 0.3).
Risque:
Le mode donne un risque plus élevé (3.475 > 2.665) pour la perte quadratique.
Exercice 15: Prédicteur de Bayes optimal (perte 0-1) ★★
Pour la classification avec perte 0-1, le prédicteur de Bayes optimal est le mode conditionnel (la classe la plus probable).
Considérons un problème de classification à 3 classes avec les probabilités conditionnelles suivantes pour un donné:
Quel est le prédicteur de Bayes optimal ?
Calculez le risque de Bayes .
Supposons qu’un classificateur prédit la classe 0 pour cet . Quel est son risque?
Situation asymétrique: Supposons que se tromper sur la classe 1 (maladie) coûte 10 fois plus cher que les autres erreurs. Définissez une matrice de coût et trouvez la prédiction optimale.
Montrez que pour la perte 0-1, aucun classificateur ne peut avoir un risque inférieur au risque de Bayes.
Solution Exercice 15
Prédicteur optimal:
Pour la perte 0-1, .
Ici, la classe 1 a la probabilité maximale (0.45), donc .
Risque de Bayes:
Même le meilleur classificateur possible se trompe 55% du temps pour ce .
Risque du classificateur sous-optimal:
Si on prédit la classe 0:
Ce classificateur a un risque plus élevé (0.75 > 0.55).
Coûts asymétriques:
Matrice de coût = coût de prédire quand la vraie classe est :
0 10 1 1 0 1 1 10 0 Risque espéré pour chaque prédiction:
Prédire 0:
Prédire 1:
Prédire 2:
Prédiction optimale: classe 1 (coût minimal 0.55).
Optimalité:
Le risque est .
Pour minimiser cette quantité, il faut maximiser , donc choisir .
Tout autre choix donne un risque plus élevé. Le prédicteur de Bayes est donc optimal par construction.
Exercice 16: SVD et facteurs de rétrécissement ★★★
La décomposition en valeurs singulières (SVD) de révèle pourquoi Ridge “rétrécit” les coefficients de manière différenciée.
Pour la matrice de données suivante, calculez la SVD :
Vérifiez que (les colonnes de sont les vecteurs propres de ).
Pour , calculez les facteurs de rétrécissement pour chaque direction .
Expliquez pourquoi la direction avec la plus petite valeur singulière est plus fortement rétrécée.
Tracez les facteurs de rétrécissement et en fonction de pour .
Solution Exercice 16
SVD:
X = np.array([[2, 1], [2, 2], [2, 3]]) U, d, Vt = np.linalg.svd(X, full_matrices=False) V = Vt.T D = np.diag(d)Résultat (approximatif):
,
,
Vérification:
XtX = X.T @ X VD2Vt = V @ D**2 @ V.T np.allclose(XtX, VD2Vt) # TrueOn peut aussi vérifier que les valeurs propres de sont et .
Facteurs de rétrécissement pour λ = 1:
Explication:
La direction 2 a une petite valeur singulière (), ce qui signifie que les données varient peu dans cette direction. L’information est donc “faible” et potentiellement bruitée. Ridge pénalise plus fortement cette direction ( vs ) pour éviter d’ajuster le bruit.
En termes de conditionnement: le rapport indique que la matrice est mal conditionnée. Ridge améliore ce conditionnement en réduisant l’effet des petites valeurs singulières.
Visualisation:
lambdas = np.linspace(0, 10, 100) s1 = d[0]**2 / (d[0]**2 + lambdas) s2 = d[1]**2 / (d[1]**2 + lambdas) plt.plot(lambdas, s1, label=f's1 (d1={d[0]:.2f})') plt.plot(lambdas, s2, label=f's2 (d2={d[1]:.2f})') plt.xlabel('λ') plt.ylabel('Facteur de rétrécissement') plt.legend()On observe que reste proche de 1 même pour modéré, tandis que décroît rapidement. C’est le rétrécissement différencié de Ridge.
Exercice 17: Conditionnement et stabilité numérique ★★★
Le nombre de conditionnement mesure la sensibilité de la solution d’un système linéaire aux perturbations.
Pour une matrice symétrique définie positive, où sont les valeurs propres.
Calculez le nombre de conditionnement de pour:
Résolvez le système pour .
Perturbez légèrement en et résolvez à nouveau. Comment change la solution?
Montrez que pour Ridge, .
Calculez le conditionnement de la matrice Ridge pour . Comparez avec le cas MCO.
Solution Exercice 17
Conditionnement de X’X:
X = np.array([[1, 1], [1, 1.001], [1, 0.999]]) A = X.T @ X eigvals = np.linalg.eigvalsh(A) kappa = eigvals.max() / eigvals.min() print(f"Conditionnement: {kappa:.0f}")Le conditionnement est très élevé (de l’ordre de 106) car les colonnes sont presque colinéaires.
Solution MCO:
y = np.array([1, 2, 3]) theta = np.linalg.solve(A, X.T @ y)La solution peut être numériquement instable.
Perturbation:
y_perturb = np.array([1.01, 2, 3]) # 1% de perturbation sur y[0] theta_perturb = np.linalg.solve(A, X.T @ y_perturb) print(f"Changement: {np.linalg.norm(theta_perturb - theta)}")Une perturbation de 1% sur peut causer un changement de plusieurs centaines de % sur . C’est le signe d’un système mal conditionné.
Conditionnement Ridge:
La matrice Ridge est . Ses valeurs propres sont (où sont les valeurs propres de ).
Pour , le numérateur et le dénominateur sont tous deux augmentés, mais le dénominateur relativement plus (puisque est petit). Le conditionnement diminue.
Comparaison:
lambda_reg = 0.1 A_ridge = A + lambda_reg * np.eye(2) eigvals_ridge = np.linalg.eigvalsh(A_ridge) kappa_ridge = eigvals_ridge.max() / eigvals_ridge.min() print(f"Conditionnement MCO: {kappa:.0f}") print(f"Conditionnement Ridge: {kappa_ridge:.0f}")Le conditionnement Ridge est beaucoup plus faible (quelques dizaines au lieu de millions), ce qui rend le système numériquement stable.
Exercice 18: Inférence bayésienne complète ★★★
L’inférence bayésienne complète calcule la distribution a posteriori des paramètres, pas seulement son mode (MAP).
Considérons un modèle de régression linéaire bayésien avec:
Vraisemblance:
Prior:
Supposons et .
Pour une seule observation , calculez la distribution a posteriori . Utilisez le fait que le produit de deux gaussiennes est une gaussienne.
Quelle est la moyenne a posteriori et la variance a posteriori ?
Calculez l’estimateur MAP et comparez avec .
Ajoutez une deuxième observation . Mettez à jour la distribution a posteriori.
Tracez les distributions a priori, a posteriori après 1 observation, et a posteriori après 2 observations. Que remarquez-vous sur l’évolution de l’incertitude?
Solution Exercice 18
Calcul de la distribution a posteriori:
La distribution a posteriori est proportionnelle à:
En prenant le log:
C’est une forme quadratique en , donc la distribution a posteriori est gaussienne.
Paramètres de la distribution a posteriori:
Pour le modèle conjugué gaussien-gaussien:
Donc .
Estimateur MAP:
Pour une distribution a posteriori gaussienne, le mode = la moyenne:
Pour un modèle gaussien conjugué, MAP = moyenne a posteriori.
Mise à jour séquentielle:
On utilise la distribution a posteriori après la première observation comme nouvel a priori:
Visualisation:
theta = np.linspace(-2, 3, 200) prior = scipy.stats.norm.pdf(theta, 0, 1) post1 = scipy.stats.norm.pdf(theta, 1.2, np.sqrt(0.2)) post2 = scipy.stats.norm.pdf(theta, 1.083, np.sqrt(1/6)) plt.plot(theta, prior, label='A priori') plt.plot(theta, post1, label='A posteriori (1 obs)') plt.plot(theta, post2, label='A posteriori (2 obs)') plt.legend()Observations:
L’a priori est large (grande incertitude)
Après 1 observation, la distribution a posteriori se concentre autour de 1.2
Après 2 observations, la distribution a posteriori se concentre davantage (variance diminue)
La moyenne a posteriori est une moyenne pondérée de l’a priori et des données
Exercice 19: Entropie croisée et divergence de Kullback-Leibler ★★★
La divergence de Kullback-Leibler (KL) mesure la différence entre deux distributions. L’entropie croisée est liée à la KL divergence.
Définitions pour des distributions discrètes (vraie) et (modèle):
Montrez que .
Pour la vraie distribution et deux modèles et , calculez et . Quel modèle est “meilleur”?
Montrez que minimiser l’entropie croisée par rapport à revient à minimiser .
Expliquez pourquoi la KL divergence n’est pas symétrique: . Calculez les deux pour l’exemple de la question 2.
Reliez ceci au maximum de vraisemblance: si est la distribution empirique des données et est le modèle, montrez que minimiser la NLV revient à minimiser .
Solution Exercice 19
Relation entropie croisée et KL:
Calcul des KL divergences:
est meilleur car .
Minimisation:
Puisque et que ne dépend pas de :
Asymétrie de la KL:
.
L’asymétrie vient du fait que la KL pénalise différemment selon quelle distribution est au numérateur du log.
Lien avec le maximum de vraisemblance:
La distribution empirique est .
Puisque est constant:
Le maximum de vraisemblance cherche le modèle le plus proche (au sens KL) de la distribution empirique.
- Murphy, K. P. (2022). Probabilistic Machine Learning: An Introduction. MIT Press. https://probml.github.io/pml-book/book1.html